Interpretation and design of hydraulic fracturing treatments

ABSTRACT

Solutions for the propagation of a hydraulic fracture in a permeable elastic rock and driven by injection of a Newtonian fluid. Through scaling, the dependence of the solution on the problem parameters is reduced to a small number of dimensionless parameters.

RELATED APPLICATION

[0001] This application claims the benefit of priority under 35 U.S.C.119(e) to U.S. Provisional Patent Application Serial No. 60/353,413,filed Feb. 1, 2002, which is incorporated herein by reference.

FIELD

[0002] The present invention relates generally to fluid flow, and morespecifically to fluid flow in hydraulic fracturing operations.

BACKGROUND

[0003] A particular class of fractures in the Earth develops as a resultof internal pressurization by a viscous fluid. These fractures areeither man-made hydraulic fractures created by injecting a viscous fluidfrom a borehole, or natural fractures such as kilometers-long volcanicdikes driven by magma coming from the upper mantle beneath the Earth'scrust. Man-made hydraulic fracturing “treatments” have been performedfor many decades, and for many purposes, including the recovery of oiland gas from underground hydrocarbon reservoirs.

[0004] Despite the decades-long practice of hydraulic fracturing, manyquestions remain with respect to the dynamics of the process. Questionssuch as: how is the fracture evolving in shape and size; how is thefracturing pressure varying with time; what is the process dependence onthe properties of the rock, on the in situ stresses, on the propertiesof both the fracturing fluid and the pore fluid, and on the boundaryconditions? Some of the difficulties of answering these questionsoriginate from the non-linear nature of the equation governing the flowof fluid in the fracture, the non-local character of the elasticresponse of the fracture, and the time-dependence of the equationgoverning the exchange of fluid between the fracture and the rock.Non-locality, non-linearity, and history-dependence conspire to yield acomplex solution structure that involves coupled processes at multiplesmall scales near the tip of the fracture.

[0005] Early modeling efforts focused on analytical solutions forfluid-driven fractures of simple geometry, either straight in-planestrain or penny-shaped. They were mainly motivated by the problem ofdesigning hydraulic fracturing treatments. These solutions weretypically constructed, however, with strong ad hoc assumptions notclearly supported by relevant physical arguments. In recent years, thelimitations of these solutions have shifted the focus of research in thepetroleum industry towards the development of numerical algorithms tomodel the three-dimensional propagation of hydraulic fractures inlayered strata characterized by different mechanical properties and/orin-situ stresses. Devising a method that can robustly and accuratelysolve the set of coupled non-linear history-dependentintegro-differential equations governing this problem will advance theability to predict and interactively control the dynamic behavior ofhydraulic fracture propagation.

[0006] For the reasons stated above, and for other reasons stated belowwhich will become apparent to those skilled in the art upon reading andunderstanding the present specification, there is a need in the art foralternate methods for modeling various behaviors of hydraulic fracturingoperations.

BRIEF DESCRIPTION OF THE DRAWINGS

[0007]FIG. 1 shows a view of a radial fluid-driven fracture with anexaggerated aperture;

[0008]FIG. 2 shows a tip of a fluid-driven fracture with lag;

[0009]FIG. 3 shows a rectangular parametric space;

[0010]FIG. 4 shows a pyramid-shaped parametric space;

[0011]FIG. 5 shows a triangular parametric space;

[0012]FIG. 6 shows a semi-infinite fluid-driven crack propagating inelastic, permeable rock;

[0013]FIG. 7 shows another triangular parametric space;

[0014]FIG. 8 shows a plane strain hydraulic fracture;

[0015]FIG. 9 shows another rectangular parametric space;

[0016]FIG. 10 shows a triangular parametric space with two trajectories;

[0017]FIG. 11 shows a graph illustrating the dependence of adimensionless fracture radius on a dimensionless toughness; and

[0018]FIG. 12 shows another triangular parametric space with twotrajectories.

DESCRIPTION OF EMBODIMENTS

[0019] In the following detailed description, reference is made to theaccompanying drawings that show, by way of illustration, specificembodiments in which the invention may be practiced. These embodimentsare described in sufficient detail to enable those skilled in the art topractice the invention. It is to be understood that the variousembodiments of the invention, although different, are not necessarilymutually exclusive. For example, a particular feature, structure, orcharacteristic described herein in connection with one embodiment may beimplemented within other embodiments without departing from the spiritand scope of the invention. In addition, it is to be understood that thelocation or arrangement of individual elements within each disclosedembodiment may be modified without departing from the spirit and scopeof the invention. The following detailed description is, therefore, notto be taken in a limiting sense, and the scope of the present inventionis defined only by the appended claims, appropriately interpreted, alongwith the full range of equivalents to which the claims are entitled. Inthe drawings, like numerals refer to the same or similar functionalitythroughout the several views.

[0020] The processes associated with hydraulic fracturing includeinjecting a viscous fluid into a well under high pressure to initiateand propagate a fracture. The design of a treatment relies on theability to predict the opening and the size of the fracture as well asthe pressure of the fracturing fluid, as a function of the properties ofthe rock and the fluid. However, in view of the great uncertainty in thein-situ conditions, it is helpful to identify key dimensionlessparameters and to understand the dependence of the hydraulic fracturingprocess on these parameters. In that respect, the availability ofsolutions for idealized situations can be very valuable. For example,idealized situations such as penny-shaped (or “radial”) fluid-drivenfractures and plane strain (often referred to as “KGD,” an acronym fromthe names of researchers) fluid-driven fractures offer promise.Furthermore, the two types of simple geometries (radial and planar) arefundamentally related to the two basic types of boundary conditionscorresponding to the fluid “point”-source and the fluid “line”-source,respectively.

[0021] Various embodiments of the present invention create opportunitiesfor significant improvement in the design of hydraulic fracturingtreatments in petroleum industry. For example, numerical algorithms usedfor simulation of actual hydraulic fracturing treatments in varyingstress environment in inhomogeneous rock mass, can be significantlyimproved by embedding the correct evolving structure of the tip solutionas described herein. Also for example, various solutions of a radialfracture in homogeneous rock and constant in-situ stress presentnon-trivial benchmark problems for the numerical codes for realistichydraulic fractures in layered rocks and changing stress environment.Also, mapping of the solution in a reduced dimensionless parametricspace opens an opportunity for a rigorous solution of an inverse problemof identification of the parameters which characterize the reservoirrock and the in-situ state of stress from the data collected duringhydraulic fracturing treatment.

[0022] Various applications of man-made hydraulic fractures includesequestration of CO₂ in deep geological layers, stimulation ofgeothermal reservoirs and hydrocarbon reservoirs, cuttings reinjection,preconditioning of a rock mass in mining operations, progressive closureof a mine roof, and determination of in-situ stresses at great depth.Injection of fluid under pressure into fracture systems at depth canalso be used to trigger earthquakes, and holds promise as a technique tocontrol energy release along active fault systems.

[0023] Mathematical models of hydraulic fractures propagating inpermeable rocks should account for the primary physical mechanismsinvolved, namely, deformation of the rock, fracturing or creation of newsurfaces in the rock, flow of viscous fluid in the fracture, andleak-off of the fracturing fluid into the permeable rock. The parametersquantifying these processes correspond to the Young's modulus E andPoisson's ratio v, the rock toughness K_(lc), the fracturing fluidviscosity μ (assuming a Newtonian fluid), and the leak-off coefficientC_(l), respectively. There is also the issue of the fluid lag λ, thedistance between the front of the fracturing fluid and the crack edge,which brings into the formulation the magnitude of far-field stressσ_(o) (perpendicular to the fracture plane) and the virgin pore pressurep_(o).

[0024] Multiple embodiments of the present invention are described inthis disclosure. Some embodiments deal with radial hydraulic fractures,and some other embodiments deal with plane strain (KGD) fractures, andstill other embodiments are general to all types of fractures. Further,different embodiments employ various scalings and various parametricspaces. For purposes of illustration, and not by way of limitation, theremainder of this disclosure is organized by different types ofparametric spaces, and various other organizational breakdowns areprovided within the discussion of the different types of parametricspaces.

[0025] I. Embodiments Utilizing a First Parametric Space

[0026] A. Radial Fractures

[0027] The problem of a radial hydraulic fracture driven by injecting aviscous fluid from a “point”-source, at a constant volumetric rate Q_(o)is schematically shown in FIGS. 1 and 2. Under conditions where the lagis negligible (λ/R<<1), determining the solution of this problemconsists of finding the aperture w of the fracture, and the net pressurep (the difference between the fluid pressure p_(f) and the far-fieldstress σ_(o)) as a function of both the radial coordinate r and time t,as well as the evolution of the fracture radius R(t). The functionsR(t), w(r,t), and p(r,t) depend on the injection rate Q_(o) and on the 4material parameters E′, μ′, K′, and C′ respectively defined as$\begin{matrix}\begin{matrix}{E^{\prime} = \frac{E}{1 - v^{2}}} & {\mu^{\prime} = {12\mu}} & {K^{\prime} = {4\left( \frac{2}{\pi} \right)^{1/2}K_{lc}}} & {C^{\prime} = {2C_{l}}}\end{matrix} & (1)\end{matrix}$

[0028] The three functions R(t), w(r,t), and p(r,t) are determined bysolving a set of equations which can be summarized as follows.

[0029] Elasticity Equation: $\begin{matrix}{w = {\frac{R}{E^{\prime}}{\int_{0}^{1}{{G\left( {{r/R},s} \right)}{p\left( {{sR},t} \right)}s\quad {s}}}}} & (2)\end{matrix}$

[0030] where G is a known elastic kernel. This singular integralequation expresses the non-local dependence of the fracture width w onthe net pressure p.

[0031] Lubrication Equation: $\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{1}{r}\frac{\partial\quad}{\partial r}\left( {{rw}^{3}\frac{\partial p}{\partial r}} \right)}} & (3)\end{matrix}$

[0032] This non-linear differential equation governs the flow of viscousincompressible fluid inside the fracture. The function g(r,t) denotesthe rate of fluid leak-off, which evolves according to $\begin{matrix}{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(r)}}}} & (4)\end{matrix}$

[0033] where t_(o) (r) is the exposure time of point r (i.e., the timeat which the fracture front was at a distance r from the injectionpoint). The leak-off law (4) is an approximation with the constant C′lumping various small scale processes (such as displacement of the porefluid by the fracturing fluid). In general, (4) can be defended underconditions where the leak-off diffusion length is small compared to thefracture length.

[0034] Global Volume Balance: $\begin{matrix}{{Q_{o}t} = {{2\pi {\int_{0}^{R}{{wr}\quad {r}}}} + {2\pi {\int_{0}^{r}{r\quad {\int_{0}^{R{(\tau)}}\quad {{g\left( {r,\tau} \right)}\quad {r}\quad {\tau}}}}}}}} & (5)\end{matrix}$

[0035] This equation expresses that the total volume of fluid injectedis equal to the sum of the fracture volume and the volume of fluid lostin the rock surrounding the fracture.

[0036] Propagation Criterion: $\begin{matrix}\begin{matrix}{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{R - r}}},} & {1 - {\frac{r}{R}{1}}}\end{matrix} & (6)\end{matrix}$

[0037] Within the framework of linear elastic fracture mechanics, thisequation embodies the fact that the fracture is always propagating andthat energy is dissipated continuously in the creation of new surfacesin rock (at a constant rate per unit surface). Note that (6) impliesthat w=0 at the tip.

[0038] Tip Conditions: $\begin{matrix}\begin{matrix}{{{w^{3}\frac{\partial p}{\partial r}} = 0},} & {r = R}\end{matrix} & (7)\end{matrix}$

[0039] This zero fluid flow rate condition (q=0) at the fracture tip isapplicable only if the fluid is completely filling the fracture(including the tip region) or if the lag is negligible at the scale ofthe fracture. Otherwise, the equations have to be altered to account forthe phenomena taking place in the lag zone as discussed below.Furthermore, the lag size λ(t) is unknown, see FIG. 2.

[0040] The formulated model for the radial fracture or similar model fora planar fracture gives a rigorous account for various physicalmechanisms governing the propagation of hydraulic fractures, however, isbased on number of assumptions which may not hold for some specificclasses of fractures. Particularly, the effect of fracturing fluidbuoyancy (the difference between the density of fracturing fluid and thedensity of the host rock) is one of the driving mechanisms of verticalmagma dykes (though, inconsequential for the horizontal disk shapedmagma fractures) is not considered in this proposal. Other processeswhich could be relevant for the hydraulic fracture propagation undercertain limited conditions which are not discussed here include aprocess zone near the fracture tip, fracturing fluid cooling andsolidification effects (as relevant to magma-driven fractures),capillarity effects at the fluid front in the fracture, and deviationsfrom the one-dimensional leak-off law.

[0041] 1. Propagation Regimes of Finite Fractures

[0042] Scaling laws for finite radial fracture driven by fluid injectedat a constant rate are considered next. Similar scaling can be developedfor other geometries and boundary conditions. Regimes with negligiblefluid lag are differentiated from regimes with non-negligible fluid lag.

[0043] a. Regimes with Negligible Fluid Lag.

[0044] Propagation of a hydraulic fracture with zero lag is governed bytwo competing dissipative processes associated with fluid viscosity andsolid toughness, respectively, and two competing components of the fluidbalance associated with fluid storage in the fracture and fluid storagein the surrounding rock (leak-off). Consequently, limiting regimes ofpropagation of a fracture can be associated with dominance of one of thetwo dissipative processes and/or dominance of one of the two fluidstorage mechanisms. Thus, four primary asymptotic regimes of hydraulicfracture propagation with zero lag can be identified where one of thetwo dissipative mechanisms and one of the two fluid storage componentsare vanishing: storage-viscosity (M), storage-toughness (K),leak-off-viscosity ({tilde over (M)}), and leak-off-toughness ({tildeover (K)}) dominated regimes. For example, fluid leak-off is negligiblecompared to the fluid storage in the fracture and the energy dissipatedin the flow of viscous fluid in the fracture is negligible compared tothe energy expended in fracturing the rock in thestorage-viscosity-dominated regime (M). The solution in thestorage-viscosity-dominated regime is given by the zero-toughness,zero-leak-off solution (K′=C′=0). As used herein, the letters M (forviscosity) and K (for toughness) are used to identify which dissipativeprocess is dominant and the symbol tilde (˜) (for leak-off) and no-tilde(for storage in the fracture) are used to identify which fluid balancemechanism is dominant.

[0045] Consider general scaling of the finite fracture which hinges ondefining the dimensionless crack opening Ω, net pressure Π, and fractureradius γ as:

w=εΩ(ρ;P ₁ ,P ₂), p=εE′Π(ρ;P ₁ ,P ₂), R=γ(P ₁ ,P ₂)L  (8)

[0046] These definitions introduce a scaled coordinate ρ=r/R(t) (0≦ρ≦1),a small number ε(t), a length scale L(t) of the same order of magnitudeas the fracture length R(t), and two dimensionless evolution parametersP₁(t) and P2(t), which depend monotonically on t. The form of thescaling (8) can be motivated from elementary elasticity considerations,by noting that the average aperture scaled by the fracture radius is ofthe same order as the average net pressure scaled by the elasticmodulus.

[0047] Four different scalings can be defined to emphasize abovedifferent primary limiting cases. These scalings yield power lawdependence of L, ε, P₁, and P₂ on time t; i.e. L˜t^(α), ε˜t^(δ), P₁˜t⁶²^(₁) , P₂˜t^(β) ^(₂) , see Table 1 for the case of a radial fracture.Furthermore, the evolution parameters can take either the meaning of atoughness (K_(m), K_({tilde over (m)})), or a viscosity (M_(k),M_({tilde over (k)})), or a storage (S_({tilde over (m)}),S_({tilde over (k)})) or a leak-off coefficient (C_(m), C_(k)). Scalingε L P₁ P₂ storage/ viscosity(M)$\left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}$

$\left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/9}$

$K_{m} = {K^{\prime}\left( \frac{t^{2}}{\mu^{\prime 5}E^{\prime 13}Q_{o}^{3}} \right)}^{1/18}$

$C_{m} = {C^{\prime}\left( \frac{{E^{\prime}}^{3}t^{7}}{\mu^{\prime 4}Q_{o}^{6}} \right)}^{1/18}$

storage/ toughness(K)$\left( \frac{K^{\prime 6}}{E^{\prime 6}Q_{o}t} \right)^{1/5}$

$\left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/5}$

$M_{k} = {\mu^{\prime}\left( \frac{Q_{o}^{3}E^{\prime 13}}{K^{\prime 18}t^{2}} \right)}^{1/5}$

$C_{k} = {C^{\prime}\left( \frac{E^{\prime 8}t^{3}}{K^{\prime 8}Q_{o}^{2}} \right)}^{1/10}$

leak-off/ viscosity(M)$\left( \frac{\mu^{\prime 4}C^{\prime 6}}{E^{\prime 4}Q_{o}^{2}t^{3}} \right)^{\frac{1}{16}}$

$\left( \frac{Q_{o}^{2}t}{C^{\prime 2}} \right)^{1/4}$

$K_{\overset{\sim}{m}} = {K^{\prime}\left( \frac{t}{E^{\prime 12}\mu^{\prime 4}C^{\prime 2}Q_{o}^{2}} \right)}^{\frac{1}{16}}$

$S_{\overset{\sim}{m}} = \left( \frac{\mu^{\prime 4}Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}t^{7}} \right)^{1/16}$

leak-off/ toughness(K)$\left( \frac{K^{\prime 8}C^{\prime 2}}{E^{\prime 8}Q_{o}^{2}t} \right)^{1/8}$

$\left( \frac{Q_{o}^{2}t}{C^{\prime 2}} \right)^{1/4}$

$M_{\overset{\sim}{k}} = {\mu^{\prime}\left( \frac{E^{\prime 12}C^{\prime 2}Q_{o}^{2}}{K^{\prime 16}t} \right)}^{1/4}$

$S_{\overset{\sim}{k}} = \left( \frac{K^{\prime 8}Q_{o}^{2}}{E^{\prime 8}C^{\prime 10}t^{3}} \right)^{1/8}$

[0048] Table 1. Small parameter ε, lengthscale L, and parameters P₁ andP₂ for the two storage scalings (viscosity and toughness) and the twoleak-off scalings (viscosity and toughness).

[0049] The regimes of solutions can be conceptualized in a rectangularparametric space MK{tilde over (K)}{tilde over (M)} shown in FIG. 3.Each of the four primary regimes (M, K, {tilde over (M)}, and {tildeover (K)}) of hydraulic fracture propagation corresponding to thevertices of the diagram is dominated by only one component of fluidglobal balance while the other can be neglected (i.e. respective P₁=0,see Table 1) and only one dissipative process while the other can beneglected (i.e. respective P₂=0, see Table 1). The solution for each ofthe primary regimes has the property that it evolves with time taccording to a power law. In particular, the fracture radius R evolvesin these regimes according to l˜t^(α) where the exponent α depends onthe regime of propagation: α=4/9,2/5,1/4,1/4 in the M-, K-, {tilde over(M)}-, {tilde over (K)}-regime, respectively. As follows from thestationary tip solution (see below), the behavior of the solution at thetip also depends on the regime of solution: Ω˜(1−ρ)^(2/3) at theM-vertex, Ω˜(1−ρ)^(5/8) at the {tilde over (M)}-vertex, andΩ˜(1−ρ)^(1/2) at the K- and {tilde over (K)}-vertices.

[0050] The edges of the rectangular phase diagram MK{tilde over(K)}{tilde over (M)} can be identified with the four secondary limitingregimes corresponding to either the dominance of one of the two fluidglobal balance mechanisms or the dominance of one of the two energydissipation mechanisms: storage-edge (MK, C_(m)=C_(k)=0), leak-off-edge({tilde over (M)}{tilde over (K)},S_({tilde over (m)})=S_({tilde over (k)})=0), viscosity-edge (M{tildeover (M)}, K_(m)=K_({tilde over (m)})=0), and K{tilde over(K)}-toughness-edge (M_(k)=M_({tilde over (k)})=0).

[0051] The regime of propagation evolves with time, since the parametersM's, K's, C's and S's depend on t. With respect to the evolution of thesolution in time, it is useful to locate the position of the state pointin the MK{tilde over (K)}{tilde over (M)} space in terms of thedimensionless times τ_(mk)=t/t_(mk),τ_({tilde over (m)}{tilde over (k)})=t/t_({tilde over (m)}{tilde over (k)}),τ_({tilde over (m)}{tilde over (m)})=t/t_({tilde over (m)}{tilde over (m)}),and τ_(k{overscore (k)})=t/t_(k{overscore (k)}) where the time scalesare defined as $\begin{matrix}{{t_{mk} = \left( \frac{\mu^{\prime 5}E^{\prime 13}Q_{o}^{3}}{K^{\prime 18}} \right)^{1/2}},} & \left( {9a} \right) \\{t_{\overset{\sim}{m}\quad \overset{\sim}{k}} = \frac{\mu^{\prime 4}E^{\prime 2}C^{\prime 2}Q_{o}^{2}}{K^{\prime 6}}} & \left( {9b} \right) \\{t_{m\quad \overset{\sim}{m}} = \left( \frac{\mu^{\prime 4}Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}} \right)^{1/7}} & \left( {9c} \right) \\{t_{k\quad \overset{\sim}{k}} = \left( \frac{\mu^{\prime 8}Q_{o}^{2}}{E^{\prime 8}C^{\prime 10}} \right)^{1/3}} & \left( {9d} \right)\end{matrix}$

[0052] Only two of these times are independent, however, sincet_({tilde over (m)}{tilde over (k)})=t_(mk) ^(8/5)t_(k{tilde over (k)})^(−3/5) and t_(m{tilde over (m)})=t_(mk) ^(8/35)t_(k{tilde over (k)})^(27/35). Note that the parameters M's, K's, C's and S's can be simplyexpressed in terms of these times according to $\begin{matrix}{{K_{m} = {M_{k}^{{- 5}/18} = \tau_{mk}^{1/9}}},{K_{\overset{\sim}{m}} = {M_{\overset{\sim}{k}}^{{- 1}/4} = \tau_{\overset{\sim}{m}\quad \overset{\sim}{k}}^{1/16}}},{C_{m} = {S_{\overset{\sim}{m}}^{{- 8}/9} = \tau_{m\quad \overset{\sim}{m}}^{7/18}}},{C_{k} = {S_{\overset{\sim}{k}}^{{- 4}/5} = \tau_{k\quad \overset{\sim}{k}}^{3/10}}}} & (10)\end{matrix}$

[0053] The dimensionless times τ's define evolution of the solutionalong the respective edges of the rectangular space MK{tilde over(K)}{tilde over (M)}. A point in the parametric space MK{tilde over(K)}{tilde over (M)} is thus completely defined by any pair combinationof these four times, say (τ_(mk), τ_(k{tilde over (k)})). The position(τ_(mk), τ_(k{tilde over (k)})) of the state point can in fact beconceptualized at the intersection of two rays, perpendicular to thestorage- and toughness-edges respectively. Furthermore, the evolution ofthe solution regime in the MK{tilde over (K)}{tilde over (M)} spacetakes place along a trajectory corresponding to a constant value of theparameter η, which is related to the ratios of characteristic times$\begin{matrix}{{\eta = \frac{E^{{\prime 11}/2}\mu^{{\prime 3}/2}C^{\prime 2}Q_{o}^{1/2}}{K^{\prime 7}}},{{{with}\frac{t_{m\quad \overset{\sim}{m}}}{t_{mk}}} = \eta},{\frac{t_{m\quad \overset{\sim}{m}}}{t_{mk}} = \eta},{\frac{t_{k\quad \overset{\sim}{k}}}{t_{mk}} = \eta^{{- 5}/3}},{\frac{t_{m\quad \overset{\sim}{m}}}{t_{k\quad \overset{\sim}{k}}} = \eta^{8/21}}} & (11)\end{matrix}$

[0054] (One of such trajectories is shown at 310 in FIG. 3).

[0055] In view of the dependence of the parameters M's, K's, C's, andS's on time, (10), the M-vertex corresponds to the origin of time, andthe {tilde over (K)}-vertex to the end of time (except for animpermeable rock). Thus, given all the problem parameters whichcompletely define the number η, the system evolves with time (say timeτ_(mk)) along a η-trajectory, starting from the M-vertex (K_(m)=0,C_(m)=0) and ending at the {tilde over (K)}-vertex(M_({tilde over (k)})=0, S_({tilde over (k)})=0). If η=0, twopossibilities exist: either the rock is impermeable (C′=0) and thesystem evolves along the storage edge from M to K, or the fluid isinviscid (μ′=0) and the system then evolves along the toughness-edgefrom K to {tilde over (K)}. If η=∞, then either K′=0 (corresponding to apre-existing discontinuity), and the system evolves along theviscosity-edge from M to {tilde over (M)}; or C′=∞ (corresponding tozero fluid storage in the fracture) and the system evolves along theleak-off-edge from the {tilde over (M)} to the {tilde over (K)}. Thuswhen η is decreasing (which can be interpreted for example as andecreasing ratio t_(m{tilde over (m)})/t_(mk)), the trajectory isattracted by the K-vertex, and when η is increasing the trajectory isattracted to the {tilde over (M)}-vertex. The dependence of the scaledsolution F can thus be expressed in the form F(ρ,τ;η), where τ is one ofthe dimensionless time, irrespective of the adopted scaling.

[0056] b. Regimes with Non-Negligible Fluid Lag.

[0057] Under certain conditions (e.g., when a fracture propagates alongpre-existing discontinuity K′=0 and confining stress σ_(o) is smallenough), the length of the lag between the crack tip and the fluid frontcannot be neglected with respect to the fracture size. In someembodiments of the present invention, fluid pressure in the lag zone canbe considered to be zero compared to the far-field stress σ_(o), eitherbecause the rock is impermeable or because there is cavitation of thepore fluid. Under these conditions, the presence of the lag brings σ_(o)in the problem description, through an additional evolution parameterP₃(t), which is denoted T_(m) in the M-scaling (or T_({tilde over (m)})in the {tilde over (M)}-scaling) and has the meaning of dimensionlessconfining stress. This extra parameter can be expressed in terms of anadditional dimensionless time as $\begin{matrix}{T_{m} = {{\tau_{om}^{1/3}\quad {with}\quad \tau_{om}} = {{\frac{t}{t_{om}}\quad {and}\quad t_{om}} = \frac{\mu^{\prime}E^{\prime 2}}{\sigma_{o}^{3}}}}} & (12)\end{matrix}$

[0058] Now the parametric space can be envisioned as the pyramidMK{tilde over (K)}{tilde over (M)}-OÕ, depicted in FIG. 4, with theposition of the state point identified by a triplet, e.g., (T_(m),K_(m),C_(k)) or (τ_(om),τ_(mk),τ_(k{tilde over (k)})). In accord with thediscussion of the zero lag case, OÕ-edge corresponds to theviscosity-dominated regime (K_(m)=K_({tilde over (m)})=0) undercondition of vanishing confining stress (T_(m)=T_({tilde over (m)})=0),where the endpoints, O- and Õ-vertices correspond to the limits ofstorage and leak-off-dominated cases.

[0059] The system evolves from the O-vertex towards the {tilde over(K)}-vertex following a trajectory which depends on all the parametersof the problem (410, FIG. 4). The trajectory depends on two numberswhich can be taken as η defined in (11) (independent of σ_(o)) andφ=t_(om)/t_(mk). It should be noted that the O-vertex from wherefracture evolution initiates is a singular point as (i) it correspondsto the infinitely fast initial fracture propagation (propagation of anunconfined fracture, σ_(o)=0, along preexisting discontinuity, K′=0)(ii) it corresponds to the infinite multitude of self-similar solutionsparameterized by the ray along which the solution trajectory is emergingfrom the O-vertex.

[0060] If φ<<1 and φ<<η (e.g. the confining stress σ₀ is “large”), thetrajectory follows essentially the OM-edge, and then from the M-vertexremains within the MK{tilde over (K)}{tilde over (M)}-rectangle.Furthermore, the transition from O to M takes place extremely morerapidly than the evolution from the M to the {tilde over (K)}-vertexalong a η-trajectory (or from M to the K-vertex if the rock isimpermeable). In other words, the parametric space can be reduced to theMK{tilde over (K)}{tilde over (M)}-rectangle, and the lag can thus beneglected if φ<<1 and φ<<η. Through this reduction in the dimensions ofthe parametric space, the M-vertex becomes the apparent starting pointof the evolution of a fluid-driven fracture without lag. The “penalty”for this reduction is a multiple boundary layer structure of thesolution near the M-vertex.

[0061] If the rock is impermeable (C′=0), the solution is restricted toevolve on the MKO face of the parametric space (see FIG. 5), from O to Kfollowing a (φ-trajectory 510. However, there is no additional timescale associated with the OK-edge and thus the transition OK takes place“rapidly” if φ>>1; this is a limiting case where the lag can beneglected, as the solution is always in the asymptotic K-regime.

[0062] 2. Structure of the Solution Near the Tip of PropagatingHydraulic Fracture

[0063] The nature of the solution near the tip of a propagatingfluid-driven fracture can be investigated by analyzing the problem of asemi-infinite fracture propagating at a constant speed V, see FIGS. 6and 7. In the following, a distinction is made between regimes/scaleswith negligible and non-negligible lag between the crack tip and thefluid front. Although a lag of a priori unknown length λ between thecrack tip and the fluid front must necessarily exist on a physicalground, as otherwise the fluid pressure at the tip has negativesingularity, there are circumstances where λ is small enough compared tothe relevant lengthscales that it can be neglected. (This issue issimilar to the use of the solutions of linear elastic fracture mechanicswhich yield “unphysical” stress singularity at the fracture tip). Inthese regimes/scales, the solution is characterized by a singularbehavior, with the nature of the singularity being a function of theproblem parameters and the scale of reference.

[0064] a. Regimes/Scales with Negligible Fluid Lag.

[0065] In view of the stationary nature of the considered tip problem,the fracture opening ŵ, net pressure {circumflex over (p)} and flow rate{circumflex over (q)} are only a function of the moving coordinate{circumflex over (x)}, see FIGS. 6 and 7. The system of equationsgoverning ŵ({circumflex over (x)}) and {circumflex over (p)}({circumflexover (x)}) can be written as $\begin{matrix}{{\hat{p} = {\frac{E^{\prime}}{4\pi}{\int_{0}^{\infty}{\frac{\hat{w}}{\hat{s}}\quad \frac{\hat{s}}{\hat{x} - \hat{s}}}}}},{{\frac{{\hat{w}}^{3}}{\mu^{\prime}}\frac{\hat{p}}{\hat{x}}} = {{V\quad \hat{w}} + {2C^{\prime}V^{1/2}{\hat{x}}^{1/2}}}},{{\lim\limits_{\hat{x}0}\frac{\hat{w}}{{\hat{x}}^{1/2}}} = \frac{K^{\prime}}{E^{\prime}}},{{\lim\limits_{\hat{x}0}{{\hat{w}}^{3}\frac{\hat{p}}{\hat{x}}}} = 0}} & (13)\end{matrix}$

[0066] The singular integral equation (13)_(a) derives from elasticity,while the Reynolds equation (13)_(b) is deduced from the Poiseuille({circumflex over (q)}=ŵ³/μ′d{circumflex over (p)}/d{circumflex over(x)}), continuity (V dŵ/d{circumflex over (x)}−d{circumflex over(q)}/d{circumflex over (x)}+ĝ=0), and Carter's leak-off laws(ĝ=C′{square root}{square root over (V/{circumflex over (x)})}).Equation (13)_(c) expresses the crack propagation criterion, while thezero flow rate condition at the tip, (13)_(d), arises from theassumption of zero lag.

[0067] Analogously to the considerations for the finite fracture, fourprimary limiting regimes of propagation of a semi-infinite fracture withzero lag can be identified where one of the two dissipative mechanismsand one of the two fluid storage components are vanishing:storage-viscosity (m), storage-toughness (k), leak-off-viscosity ({tildeover (m)}), and leak-off-toughness ({tilde over (k)}) dominated regimes.Each of the regimes correspond to the respective vertex of therectangular parametric space of the semi-infinite fracture. However, inthe context of the semi-infinite fracture, the storage-toughness (k) andleak-off-toughness ({tilde over (k)}) dominated regimes are identicalsince the corresponding zero viscosity (μ′=0) solution of (13) isindependent of the balance between the fluid storage and leak-off, andis given by the classical linear elastic fracture mechanics (LEFM)solution ŵ=(K′/E′){circumflex over (x)}^(1/2) and {circumflex over(p)}=0. Therefore, the toughness edge k{tilde over (k)} of therectangular parameteric space for the semi-infinite fracture collapsesinto a point, which can be identified with either k- or {tilde over(k)}-vertex, and the rectangular space itself into the triangularparametric space mk{tilde over (m)}, see FIG. 7.

[0068] The primary storage-viscosity, toughness, and leak-off-viscosityscalings associated with the three primary limiting regimes (m, k or{tilde over (k)}, and {tilde over (m)}) are as follows $\begin{matrix}{{{\hat{\xi}}_{m,k,\overset{\sim}{m}} = \frac{\hat{x}}{l_{m,k,\overset{\sim}{m}}}},{{\hat{\Omega}}_{m,k,\overset{\sim}{m}} = \frac{\hat{w}}{l_{m,k,\overset{\sim}{m}}}},{\underset{m,k,\overset{\sim}{m}}{\overset{\quad}{\hat{\prod}}}\quad {= \frac{\hat{p}}{E^{\prime}}}},{{\hat{\Psi}}_{m,k,\overset{\sim}{m}} = \frac{\hat{q}}{V\quad l_{m,k,\overset{\sim}{m}}}}} & (14)\end{matrix}$

[0069] where the three lengthscales l_(m), l_(k) andl_({tilde over (m)}) are defined as l_(m)=μ′V/E′, l_(k)=(K′/E′)²,l_({tilde over (m)}=V) ^(1/3)(2μ′C′)^(2/3)/E′^(2/3). The solutionF = {Ω̂, Π̂}

[0070] in the various scalings can be shown to be of the form{circumflex over (F)}_(m)({circumflex over (ξ)}_(m); c_(m),k_(m)),{circumflex over (F)}_(k)({circumflex over(ξ)}_(k);m_(k),m_({tilde over (k)})), {circumflex over(F)}_({tilde over (m)})(ξ_({tilde over (m)});s_({tilde over (m)}),k_({tilde over (m)})),with the letters m's, k's, s's and c's representing dimensionlessviscosity, toughness, storage, and leak-off coefficient, respectively.$\begin{matrix}{{k_{m} = {m_{k}^{{- 1}/2} = \left( {l_{k}l_{m}} \right)^{1/2}}},{k_{\overset{\sim}{m}} = {m_{\overset{\sim}{k}}^{{- 1}/3} = \left( {l_{k}l_{\overset{\sim}{m}}} \right)^{1/2}}},{c_{m} = {s_{\overset{\sim}{m}}^{{- 3}/2} = \left( {l_{\overset{\sim}{m}}l_{m}} \right)^{3/2}}}} & (15)\end{matrix}$

[0071] For example, a point in the mk{tilde over (m)} ternary diagramcorresponds to a certain pair (k_(m), c_(m)) in the viscosity scaling,with the m-vertex corresponding to c_(m)=0 and k_(m)=0. The vertexsolutions (denoted by the subscript ‘0’) are given by $\begin{matrix}\begin{matrix}\begin{matrix}{{{\hat{\Omega}}_{m0} = {\beta_{m0}{\hat{\xi}}_{m}^{2/3}}},{{\underset{m0}{\overset{\quad}{\hat{\prod}}}\quad {= {\delta_{m0}{\hat{\xi}}_{m}^{{- 1}/3}}}};}} \\{{{\hat{\Omega}}_{k0} = {\hat{\xi}}_{k}^{1/2}},{{\underset{k0}{\overset{\quad}{\hat{\prod}}}\quad {= 0}};{and}}}\end{matrix} \\{{{\hat{\Omega}}_{\overset{\sim}{m}0} = {\beta_{\overset{\sim}{m}0}\quad {\hat{\xi}}_{\overset{\sim}{m}}^{5/8}}},{\underset{\overset{\sim}{m}0}{\overset{\quad}{\hat{\prod}}}\quad {= {\delta_{\overset{\sim}{m}0}{\hat{\xi}}_{\overset{\sim}{m}}^{{- 3}/8}}}}}\end{matrix} & (16)\end{matrix}$

[0072] with β_(m0)=2^(1/3)3^(5/6), δ_(m0)=−6^(−2/3),β_({tilde over (m)}0)≅2.534, δ_({tilde over (m)}0)≅−0.164. Thus whenthere is only viscous dissipation (edge m{tilde over (m)} correspondingto fracture propagation along preexisting discontiuity K′=0) the tipbehavior is of the form ŵ˜{circumflex over (x)}^(2/3), {circumflex over(p)}˜−{circumflex over (x)}^(−1/3) in the storage-dominated case,m-vertex, (impermeable rock C′=0) and of the form ŵ˜{circumflex over(x)}^(5/8), {circumflex over (p)}˜−{circumflex over (x)}^(−3/8) in theleak-off dominated case, {tilde over (m)}-vertex. On the other hand, thek-vertex pertains to a fracture driven by an inviscid fluid (μ′=0); thisvertex is associated with the classical tip solution of linear elasticfracture mechanics ŵ˜{circumflex over (x)}^(1/2). The general case of afluid-driven fracture with no leak-off (C′=0) or negligible storagenaturally corresponds to the mk- or {tilde over (m)}k-edges,respectively. However, a more general interpretation of the mk{tildeover (m)} parametric space can be seen by expressing the numbers m's,k's, s's, and c's in terms of a dimensionless velocity v, and aparameter {circumflex over (η)} which only depends on the parameterscharacterizing the solid and the fluid $\begin{matrix}{{v = {\frac{l_{m}}{l_{k}} = \frac{V}{V_{*}}}},{\hat{\eta} = {\frac{l_{m}l_{k}^{2}}{l_{\overset{\sim}{m}}^{3}} = \frac{4E^{\prime 3}\mu^{\prime}C^{\prime 2}}{K^{\prime 4}}}}} & (17)\end{matrix}$

[0073] where V*=K′²/μ′E is a characteristic velocity. Hence,k_(m)=v^(−1/2), k_({tilde over (m)})={circumflex over(η)}^(−1/6)v^(−1/6), c_(m)={circumflex over (η)}^(1/2)v⁻¹. The aboveexpressions indicate that the solution moves from the m-vertex towardsthe k-vertex with decreasing dimensionless velocity v, along atrajectory which depends only on {circumflex over (η)}. With increasing{circumflex over (η)}, the trajectory is pulled towards the {tilde over(m)}-vertex. Since the tip velocity of a finite fracture decreases withtime (at least under constant injection rate), the tip solutioninterpreted from this stationary solution is seen to evolve with time.In other words, as the length scales l_(m) and l_({tilde over (m)})evolve with time, the nature of the solution in the tip region at agiven physical scale evolves accordingly.

[0074] The solution along the edges of the mkm-triangle, namely, theviscosity mm-edge (k_(m)=0 ork_({tilde over (m)})=0), the storagemk-edge (c_(m)=0 or m_({tilde over (k)})=0), and the {tilde over(m)}k-edge (s_({tilde over (m)})=0 or m_(k)=0) has been obtained both inthe form of series expansion in the neighborhood of the vertices andnumerically for finite values of the non-zero parameters. These resultswere obtained in part by recognizing that the solution can be furtherresealed along each edge to eliminate the remaining parameter. Forexample, the tip solution along the mk-edge, which is governed byparameter k_(m) in the m-scaling, upon rescaling to the mixed scalingcan be expressed as {circumflex over (F)}_(mk)({circumflex over(ξ)}_(mk)) where {circumflex over (ξ)}_(mk)={circumflex over (x)}/l_(mk)with l_(mk)=l_(k) ³/l_(m) ².

[0075] The m{tilde over (m)}-, mk-, and {tilde over (m)}k-solutionsobtained so far give a glimpse on the changing structure of the tipsolution at various scales, and how these scales change with the problemparameters, in particular with the tip velocity v. Consider for examplethe mk-solution (edge of the triangle corresponding to the case ofimpermeable rock) for the opening {circumflex over (Ω)}_(mk)({circumflexover (ξ)}_(mk)), with {circumflex over (Ω)}_(mk)=k_(m) ⁻⁴{circumflexover (Ω)}_(m)=m_(k){circumflex over (Ω)}_(k). Expansion of the{circumflex over (Ω)}_(mk) at {circumflex over (ξ)}_(mk)=0 and at{circumflex over (ξ)}_(mk)=∞ is of the form $\begin{matrix}{{{\hat{\Omega}}_{mk} = {{{\beta_{m0}{\hat{\xi}}_{mk}^{2/3}} + {\beta_{mk1}{\hat{\xi}}_{mk}^{h}} + {{O\left( {\hat{\xi}}_{mk}^{{- 1}/3} \right)}\quad {at}\quad {\hat{\xi}}_{mk}}} = \infty}},{{\hat{\Omega}}_{mk} = {{{\hat{\xi}}_{mk}^{1/2} + {\beta_{k\quad {m1}}{\hat{\xi}}_{mk}} + {O\left( {\hat{\xi}}_{mk}^{3/2} \right)\quad {at}\quad {\hat{\xi}}_{mk}}} = 0}}} & (18)\end{matrix}$

[0076] The exponent h≅0.139 in the “alien” term {circumflex over(ξ)}_(mk) ^(h) of the far-field expansion (18)₁ is the solution ofcertain transcendental equation obtained in connection withcorresponding boundary layer structure. In this case, the boundary layerarises because ŵ˜{circumflex over (x)}^(1/2) near {circumflex over(x)}=0 if K′>0, but ŵ˜{circumflex over (x)}^(2/3) when K′=0. Thebehavior of the mk-solution at infinity corresponds to the m-vertexsolution. The mk-solution shows that {circumflex over (Ω)}_(mk)Ω̂_(mk) ≃ β_(m0)ξ̂_(m)^(2/3)

[0077] for {circumflex over (ξ)}_(mk)>{circumflex over (ξ)}_(mk∞), with{circumflex over (ξ)}_(mk∞)=O(1), with the consequence that there willbe corresponding practical range of parameters for which the globalsolution for C′=0 is characterized by the m-vertex asymptotic behaviorŵ˜{circumflex over (x)}^(2/3), {circumflex over (p)}˜−{circumflex over(x)}^(−1/3) (viscous dissipation only), although ŵ˜{circumflex over(x)}^(1/2) in a very small region near the tip. Taking for example V≅1m/s, E≅10³ MPa, μ′≅10⁻⁶ MPa·s, K′≅1 MPa·m^(1/2), and C′=0, thenl_(mk)≅10⁻² m. Hence, at distance larger than 10⁻² m, the solutionbehaves as if the impermeable rock has no toughness and there is onlyviscous dissipation. As discussed further below, the m-vertex solutiondevelops as an intermediate asymptote at some small distance from thetip in the finite fracture, provided the lengthscale l_(mk) is muchsmaller than the fracture dimension R.

[0078] b. Regimes/Scales with Non-Negligible Fluid Lag.

[0079] The stationary problem of a semi-infinite crack propagating atconstant velocity V is now considered, taking into consideration theexistence of a lag of a priori unknown length λ between the crack tipand the fluid front, see FIG. 2. First, considerations are restricted toimpermeable rocks. In this case, the tip cavity is filled with fluidvapors, which can be assumed to be at zero pressure.

[0080] This problem benefits from different scalings in part because thefar-field stress σ_(o), directly influences the solution, through thelag. Consider for example the mixed stress/storage/viscosity scaling(om) $\begin{matrix}{{{\hat{\xi}}_{om} = {\frac{\hat{x}}{l_{om}} = {ɛ^{3}{\hat{\xi}}_{m}}}},{{\hat{\Omega}}_{om} = {ɛ^{2}{\hat{\Omega}}_{m}}},{\hat{\prod\limits_{om}}\quad {= {ɛ^{- 1}\hat{\prod\limits_{m}}}}},{{{with}\quad l_{om}} = {ɛ^{- 3}l_{m}}},{ɛ = \frac{\sigma_{o}}{E^{\prime}}}} & (19)\end{matrix}$

[0081] It can be shown that the solution is of the form {circumflex over(F)}_(om)({circumflex over (ξ)}_(om);k_(om)), where k_(om)=ε^(1/2)k_(m)is the dimensionless toughness in this new scaling; {circumflex over(F)}_(om) behaves according to the k-vertex asymptote near the tip(Ω̂_(om) ≃ k_(om)ξ̂_(om)^(1/2)  near  ξ̂_(om) = 0)

[0082] and to the m-vertex asymptote far away from the tip ({circumflexover (Ω)}_(om) Ω̂_(om) ≃ β_(om)ξ̂_(om)^(2/3)  

[0083] for {circumflex over (ξ)}_(om)>>1). The scaled lagλ_(om)=λ/l_(om) continuously decreases with k_(om), from a maximum valueλ_(mso)≅0.36 reached either when K′=0 or σ_(o)=0. The decrease of λ_(om)with k_(om) becomes exponentially fast for large toughness (practicallywhen k_(om)>4). Furthermore, analysis of the solution indicates that{circumflex over (F)}_(om)({circumflex over (ξ)}_(om);k_(om)) can berescaled into {circumflex over (F)}_(mk)({circumflex over (ξ)}_(mk)) forlarge toughness (k_(om)>4) $\begin{matrix}{{{\hat{\xi}}_{mk} = {k_{om}^{- 6}{\hat{\xi}}_{om}}},{{\hat{\Omega}}_{mk} = {k_{om}^{- 4}{\hat{\Omega}}_{om}}},{{\hat{\Pi}}_{mk} = {k_{om}^{2}{\hat{\Pi}}_{mk}}}} & (20)\end{matrix}$

[0084] These considerations show that within the context of thestationary tip solution the fluid lag becomes irrelevant at the scalesof interest if k_(om)>4, and can thus be assumed to be zero (with theimplication that {circumflex over (q)}=0 at the tip, which leads to asingularity of the fluid pressure.) Also, the solution becomesindependent of the far-field stress σ_(o) when k_(om)>4 (except as areference value of the fluid pressure) and it can be mapped within themk{tilde over (m)} parametric space introduced earlier.

[0085] In permeable rocks, pore fluid is exchanged between the tipcavity and the porous rock and flow of pore fluid within the cavity istaking place. The fluid pressure in the tip cavity is thus unknown andfurthermore not uniform. Indeed, pore fluid is drawn in by suction atthe tip of the advancing fracture, and is reinjected to the porousmedium behind the tip, near the interface between the two fluids. (Porefluid must necessarily be returning to the porous rock from the cavity,as it would otherwise cause an increase of the lag between thefracturing fluid and the tip of the fracture, and would thus eventuallycause the fracture to stop propagating). Only elements of the solutionfor this problem exists so far, in the form of a detailed analysis ofthe tip cavity under the assumption that ŵ˜{circumflex over (x)}^(1/2)in the cavity.

[0086] This analysis shows that the fluid pressure in the lag zone canbe expressed in terms of two parameters: a dimensionless fracturevelocity {overscore (v)}=Vλ/c and a dimensionless rock permeabilityζ=kE′³/(λ^(1/2)K′³), where k and c denote respectively the intrinsicrock permeability and diffusivity. Furthermore, the solution is boundedby two asymptotic regimes: drained with the fluid pressure in the lagequilibrated with the ambient pore pressure p_(o) ({overscore (v)}<<1and {overscore (ζ)}>>1), and undrained with the fluid pressurecorresponding to its instantaneous (undrained) value at the movingfracture tip $\begin{matrix}{p_{f{({tip})}} = {p_{o} - {\frac{1}{2}\frac{K^{\prime}}{E^{\prime}}\frac{\mu_{o}}{k}\sqrt{\pi \quad {cV}}}}} & (21)\end{matrix}$

[0087] where μ_(o) is the viscosity of the pore fluid. The aboveexpression for p_(f(tip)) indicates that pore fluid cavitation can takeplace in the lag. Analysis of the regimes of solution suggests that thepore fluid pressure in the lag zone drop below cavitation limit in awide range of parameters relevant for propagation of hydraulic fracturesand magma dykes, implying a net-pressure lag condition identical to theone for impermeable rock. This condition allows one to envision theparametric space for the tip problem in the general case of thepermeable rock (leak-off) and the lag (finiteness of the confiningstress) as the pyramid mk{tilde over (m)}-oõ, where similarly to thecase of the finite fracture, see FIG. 4, vertices o- and õ-correspond tothe limits of storage and leak-off dominated cases under conditions ofvanishing toughness and confining stress. The stationary tip solutionnear the om- and õ{tilde over (m)}-edges behaves as k-vertex asymptote(ŵ˜{circumflex over (x)}^(1/2)) near the tip and as the m-vertex(ŵ˜{circumflex over (x)}^(2/3)) and m-vertex (ŵ˜{circumflex over(x)}^(5/8)) asymptote, respectively, far away from the tip.

[0088] 3. Local Tip and Global Structure of the Solution

[0089] The development of the general solution corresponding to thearbitrary η-trajectory in the MK{tilde over (K)}{tilde over (M)}rectangle (or (η,φ)-trajectory in the MK{tilde over (K)}{tilde over(M)}-OÕ pyramid) is aided by understanding the asymptotic behavior ofthe solution in the vicinities of the rectangle (pyramid) vertices andedges. These asymptotic solutions can be obtained (semi-) analyticallyvia regular or singular perturbation analysis. Construction of thosesolutions to the next order in the small parameter(s) associated withthe respective edge (or vertex) can identify the physically meaningfulrange of parameters for which the fluid-driven fracture propagates inthe respective asymptotic regime (and thus can be approximated by therespective edge (vertex) asymptotic solution). Since the solutiontrajectory evolves with time from M-vertex to the {tilde over(K)}-vertex inside of the MK{tilde over (K)}{tilde over (M)}-rectangle(or generally, from the O-vertex to the {tilde over (K)}-vertex insideof the MK{tilde over (K)}{tilde over (M)}-OÕ pyramid), it is helpful tohave valid asymptotic solutions developed in the vicinities of thesevertices. The solution in the vicinity of the some of the vertices(e.g., O, K, and {tilde over (K)}) is a regular perturbation problem,which has been solved for the K-vertex along the MK- and KO-edge of thepyramid. The solution in the vicinity of the M-vertex is challengingsince it constitutes a singular perturbation problem for a system ofnon-linear, non-local equations in more than one small parameter,namely, K_(m) (along the MK-edge), C_(m) (along the MM-edge), and,generally, E_(m)=T_(m) ⁻¹ (along the MO-edge), given that the nature ofthe tip singularity changes with a small perturbation from zero in anyof these parameters. Indeed, solution at M-vertex is given by thezero-toughness (K_(m)=0), zero leak-off (C_(m)=0), zero-lag (E_(m)=T_(m)⁻¹=0) solution which near tip behavior is given by the m-vertex tipsolution, Ω_(m)˜(1−ρ)^(2/3) and Π_(m)˜−(1−ρ)^(−1/3). Small perturbationof the M-vertex in either toughness K_(m), or leak-off C_(m), or lagE_(m) changes the nature of the near tip behavior to either thetoughness asymptote Ω_(m)˜(1−ρ)^(1/2), or the leak-off asymptoteΩ_(m)˜(1−ρ)^(5/8), or the lag asymptote Ω_(m)˜(1−ρ)^(3/2), respectively.This indicates the emergence of the near tip boundary layer (BL) whichincorporates arising toughness singularity and/or leak-off singularityand/or the fluid lag. If the perturbation is small enough, there existsa lengthscale intermediate to the fracture length and the BL “thickness”where the outer solution (i.e. the solution away from the fracture tip)and the BL solution (given by the stationary tip solution discussedabove) can be matched to form the composite solution uniformly validalong the fracture. Matching requires that the asymptotic expansions ofthe outer and the BL solutions over the intermediate lengthscale areidentical.

[0090] As an illustration, the non-trivial structure of the globalsolution in the vicinity of the M-vertex along the MK-edge (i.e.,singular perturbation problem in K_(m), while C_(m)=E_(m)=0) is nowoutlined, corresponding to the case of a fracture in impermeable rockand large confining stress (or time). The outer expansion for Ω, Π, anddimensionless fracture radius γ are perturbation expansions in powers ofK_(m)^(b),

[0091] b>0. Here the matching not only gives the coefficients in theexpansion, but also determines the exponent b. It can be shown that thetip solution along the mk-edge (18) corresponds to the O(1) term in theinner (boundary layer) expansion at the tip. The inner and outer(global) scaling for the radial fracture are related as $\begin{matrix}{{{{1 - \rho} = {\frac{81}{16\quad \gamma_{m0}^{3}}K_{m}^{6}{\hat{\xi}}_{mk}}},{\Omega_{m} = {\frac{9}{4\gamma_{m0}}K_{m}^{4}{\hat{\Omega}}_{mk}}},{\Pi_{m} = {\frac{4\gamma_{m0}}{9}K_{m}^{- 2}{\hat{\Pi}}_{mk}}}}\quad} & (22)\end{matrix}$

[0092] where γ_(mo) is the O(1) term of the outer expansion for γ givenby the M-vertex solution (K_(m)=C_(m)=E_(m)=0). Using the asymptoticexpression (18) together with the scaling (22), one finds that the outerand inner solutions match under the condition K_(m)⁶1.

[0093] <<1. Then the leading order inner and outer solutions form asingle composite solution of O(1) uniformly valid along the fracture.That is, to leading order there is a lengthscale intermediate to the tipboundary layer thickness K_(m)⁶R

[0094] R and the fracture radius R, over which the inner and outersolutions posses the same intermediate asymptote, corresponding to them-vertex solution (16)₁. This solution structure corresponds to theouter zero-toughness solution valid on the lengthscale of the fracture,and thin tip boundary layer given by the mk-edge solution.

[0095] To leading order the condition K_(m)⁶1

[0096] <<1 is merely a condition for the existence of the boundary layersolution. In order to move away from the M-vertex solution away from thetip, one has to determine the exponent b in the next term in theasymptotic expansion. From this value of b we determine the asymptoticvalidity of the approximation. This can be obtained from the next-ordermatching between the near tip asymptote in the outer expansion and theaway from tip behavior of the inner solution, see (18). Here thematching to the next order of the outer and inner solutions does notrequire the next-order inner solution, as the next order outer solutionis matched with the leading order term of the inner solution. The latterappears to be a consequence of the non-local character of theperturbation problem. Then using (18) an expression for the exponentb=4−6h is obtained which yields b≅3.18 and consequently the next ordercontribution in the asymptotic expansion away from the tip. The range ofdimensionless toughness in which fracture global (outer) solution can beapproximated by the M-vertex solution is, therefore, given byK_(m)^(3.18)1.

[0097] <<1.

[0098] B. Plane Strain (KGD) Fractures

[0099] The problem of a KGD hydraulic fracture driven by injecting aviscous fluid from a “point”-source, at a constant volumetric rate Q_(o)is schematically shown in FIG. 8. Under conditions where the lag isnegligible, determining the solution of this problem consists of findingthe aperture w of the fracture, and the net pressure p (the differencebetween the fluid pressure p_(f) and the far-field stress σ_(o)) as afunction of both the coordinate x and time t, as well as the evolutionof the fracture radius l(t). The functions l(t), w(x,t), and p(x,t)depend on the injection rate Q_(o) and on the 4 material parameters E′,μ′, K′, and C′ respectively defined as $\begin{matrix}{E^{\prime} = {{\frac{E}{1 - v^{2}}\quad \mu^{\prime}} = {{12\quad \mu \quad K^{\prime}} = {{4\left( \frac{2}{\pi} \right)^{1/2}K_{Ic}\quad C^{\prime}} = {2C_{l}}}}}} & (23)\end{matrix}$

[0100] The three functions l(t), w(x,t), and p(x,t) are determined bysolving a set of equations which can be summarized as follows.

[0101] Elasticity Equation: $\begin{matrix}{{p\left( {x,t} \right)} = {{- \frac{E^{\prime}}{4\quad \pi}}{\int_{- l}^{l}{\frac{\partial{w\left( {s,t} \right)}}{\partial s}\frac{ds}{s - x}}}}} & (24)\end{matrix}$

[0102] This singular integral equation expresses the non-localdependence of the fracture width w on the net pressure p.

[0103] Lubrication Equation: $\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{\partial\quad}{\partial x}\left( {w^{3}\frac{\partial p}{\partial x}} \right)}} & (25)\end{matrix}$

[0104] This non-linear differential equation governs the flow of viscousincompressible fluid inside the fracture. The function g(x,t) denotesthe rate of fluid leak-off, which evolves according to $\begin{matrix}{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(x)}}}} & (26)\end{matrix}$

[0105] where t_(o)(x) is the exposure time of point x (i.e., the time atwhich the fracture front was at a distance x from the injection point).

[0106] Global Volume Balance: $\begin{matrix}{{Q_{o}t} = {{2{\int_{0}^{l}{w{x}}}} + {2{\int_{0}^{t}{\int_{0}^{l{(\tau)}}{{g\left( {x,\tau} \right)}{x}{\tau}}}}}}} & (27)\end{matrix}$

[0107] This equation expresses that the total volume of fluid injectedis equal to the sum of the fracture volume and the volume of fluid lostin the rock surrounding the fracture.

[0108] Propagation Criterion: $\begin{matrix}{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{l - x}}},{1 - {\frac{x}{l}{\operatorname{<<}1}}}} & (28)\end{matrix}$

[0109] Within the framework of linear elastic fracture mechanics, thisequation embodies the fact that the fracture is always propagating andthat energy is dissipated continuously in the creation of new surfacesin rock (at a constant rate per unit surface). Note that (28) impliesthat w=0 at the tip.

[0110] Tip Conditions: $\begin{matrix}{{{w^{3}\frac{\partial p}{\partial x}} = 0},{x = l}} & (29)\end{matrix}$

[0111] This zero fluid flow rate condition (q=0) at the fracture tip isapplicable only if the fluid is completely filling the fracture(including the tip region) or if the lag is negligible at the scale ofthe fracture.

[0112] 1. Propagation Regimes of a KGD Fracture

[0113] Propagation of a hydraulic fracture with zero lag is governed bytwo competing dissipative processes associated with fluid viscosity andsolid toughness, respectively, and two competing components of the fluidbalance associated with fluid storage in the fracture and fluid storagein the surrounding rock (leak-off). Consequently, the limiting regimesof propagation of a fracture can be associated with the dominance of oneof the two dissipative processes and/or the dominance of one of the twofluid storage mechanisms. Thus, four primary asymptotic regimes ofhydraulic fracture propagation with zero lag can be identified where oneof the two dissipative mechanisms and one of the two fluid storagecomponents are vanishing: storage-viscosity (M), storage-toughness (K),leak-off-viscosity ({tilde over (M)}), and leak-off-toughness ({tildeover (K)}) dominated regimes. For example, fluid leak-off is negligiblecompared to the fluid storage in the fracture and the energy dissipatedin the flow of viscous fluid in the fracture is negligible compared tothe energy expended in fracturing the rock in thestorage-viscosity-dominated regime (M). The solution in thestorage-viscosity-dominated regime is given by the zero-toughness,zero-leak-off solution (K′=C′=0).

[0114] Consider the general scaling of the finite fracture, which hingeson defining the dimensionless crack opening Ω, net pressure Π, andfracture radius γ as

w=εLΩ(ξ;P ₁ ,P ₂), p=εE′Π(ξ;P ₁ ,P ₂), l=γ(P ₁ ,P ₂)L  (30)

[0115] With these definitions, we have introduced the scaled coordinateξ=x/l(t) (0≦ξ≦1), a small number ε(t), a length scale L(t) of the sameorder of magnitude as the fracture length l(t), and two dimensionlessevolution parameters P₁(t) and P₂(t), which depend monotonically on t.The form of the scaling (30) can be motivated from elementary elasticityconsiderations, by noting that the average aperture scaled by thefracture length is of the same order as the average net pressure scaledby the elastic modulus.

[0116] Four different scalings can be defined to emphasize abovedifferent primary limiting cases. These scalings yield power lawdependence of L, ε, P₁, and P₂ on time t; i.e. L˜t^(α), ε˜t^(δ, P)₁˜t^(β) ^(₂) , P₂˜t^(β) ^(₂) , see Table 2 for the case of a radialfracture. Furthermore, the evolution parameters can take either themeaning of a toughness (K_(m), K_({tilde over (m)})), or a viscosity(M_(k), M_({tilde over (k)})), or a storage (S_({tilde over (m)}),S_({tilde over (k)})), or a leak-off coefficient (C_(m), C_(k)). TABLE 2Small parameter ε, lengthscale L, and parameters P₁ and P₂ for the twostorage scalings (viscosity and toughness) and the two leak-off scalings(viscosity and toughness). Scaling ε L P₁ P₂ storage/ viscosity(M)$\left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}$

$\left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/6}$

$K_{m} = {K^{\prime}\left( \frac{1}{E^{\prime 3}\mu^{\prime}Q_{o}} \right)}$

$C_{m} = {C^{\prime}\left( \frac{E^{\prime}t}{\mu^{\prime}Q_{o}^{3}} \right)}^{1/6}$

storage/ toughness(K)$\left( \frac{K^{\prime 4}}{E^{\prime 4}Q_{o}t} \right)^{1/3}$

$\left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/3}$

$M_{k} = {\mu^{\prime}\left( \frac{E^{\prime 3}Q_{o}}{K^{\prime 4}} \right)}$

$C_{k} = {C^{\prime}\left( \frac{E^{\prime 4}t}{K^{\prime 4}Q_{o}^{2}} \right)}^{1/6}$

leak-off/ viscosity(M)$\left( \frac{\mu^{\prime}C^{\prime 2}}{E^{\prime}Q_{o}t} \right)^{1/4}$

$\frac{Q_{o}t^{1/2}}{C^{\prime}}$

$K_{\overset{\sim}{m}} = K_{m}$

$S_{\overset{\sim}{m}} = \left( \frac{\mu^{\prime}Q_{o}^{3}}{E^{\prime}C^{\prime 6}t} \right)^{1/4}$

leak-off/ toughness(K)$\left( \frac{K^{\prime 4}C^{\prime 2}}{E^{\prime 4}Q_{o}^{2}t} \right)^{1/4}$

$\frac{Q_{o}t^{1/2}}{C^{\prime}}$

$M_{\overset{\sim}{k}} = M_{k}$

$S_{\overset{\sim}{k}} = \left( \frac{K^{\prime 4}Q_{o}^{2}}{E^{\prime 4}C^{\prime 6}t} \right)^{1/4}$

[0117] The regimes of solutions can be conceptualized in a rectangularphase diagram MK{tilde over (K)}{tilde over (M)} shown in FIG. 9. Eachof the four primary regimes (M, K, {tilde over (M)}, and {tilde over(K)}) of hydraulic fracture propagation corresponding to the vertices ofthe diagram is dominated by only one component of fluid global balancewhile the other can be neglected (i.e. respective P₁=0, see Table 1) andonly one dissipative process while the other can be neglected (i.e.respective P₂=0, see Table 1). As follows from the stationary tipsolution, the behavior of the solution at the tip also depends on theregime of solution: Ω˜(1−ρ)^(2/3) at the M-vertex, Ω˜(1−ρ)^(5/8) at the{tilde over (M)}-vertex, and Ω˜(1−ρ)^(1/2) at the K- and {tilde over(K)}-vertices.

[0118] The edges of the rectangular phase diagram MK{tilde over(K)}{tilde over (M)} can be identified with the four secondary limitingregimes corresponding to either the dominance of one of the two fluidglobal balance mechanisms or the dominance of one of the two energydissipation mechanisms: storage-edge (MK, C_(m)=C_(k)=0), leak-off-edge({tilde over (M)}{tilde over (K)},S_({tilde over (m)})=S_({tilde over (k)})=0), viscosity-edge (M{tildeover (M)}, K_(m)=K_({tilde over (m)})=0), and K{tilde over(K)}-toughness-edge (M_(k)=M_({tilde over (k)})=0). The solution alongthe storage-edge MK and along the leak-off-edge {tilde over (M)}{tildeover (K)} has the property that it evolves with time t according to apower law, i.e., according to l˜t^(α) where the exponent α depends onthe regime of propagation: α=⅔ on the storage-edge MK and α=½ onleak-off-edge {tilde over (M)}{tilde over (K)}.

[0119] The regime of propagation evolves with time from the storage-edgeto the leak-off edge since the parameters C's and S's depend on t, butnot K's and M's. With respect to the evolution of the solution in time,it is useful to locate the position of the state point in the MK{tildeover (K)}{tilde over (M)} space in terms of η which is a power of any ofthe parameters K's and M's and a dimensionless time, eitherτ_(m{tilde over (m)})=t/t_(m{tilde over (m)}) orτ_(k{tilde over (k)}=t/t) _(k{tilde over (k)}) where $\begin{matrix}{{\eta = \frac{K^{\prime 4}}{E^{\prime 3}\mu^{\prime}Q_{o}}},{t_{m\quad \overset{\sim}{m}} = \frac{\mu^{\prime}Q_{o}^{3}}{E^{\prime}C^{\prime 6}}},{t_{k\quad \overset{\sim}{k}} = \frac{K^{\prime 4}Q_{o}^{2}}{E^{\prime 4}C^{\prime 6}}}} & (31)\end{matrix}$

[0120] also noting that τ_(m{tilde over (m)})=ηt_(k{tilde over (k)})since $\begin{matrix}{\frac{t_{k\quad \overset{\sim}{k}}}{t_{m\quad \overset{\sim}{m}}} = \eta} & (32)\end{matrix}$

[0121] The parameters M's, K's, C's and S's can be expressed in terms ofη and τ_(m{tilde over (m)})(or τ_(k{tilde over (k)})) according to

K _(m) =K _({tilde over (m)})=η^(1/4) , M _(k) =M_({tilde over (k)})=η⁻¹  (33)

[0122] $\begin{matrix}{{C_{m} = \tau_{m\quad \overset{\sim}{m}}^{1/6}},{C_{k} = {{\eta^{{- 1}/6}\tau_{m\quad \overset{\sim}{m}}^{1/6}} = \tau_{k\quad \overset{\sim}{k}}^{1/6}}}} & (34)\end{matrix}$

$\begin{matrix}{{S_{\overset{\sim}{m}} \approx \tau_{m\quad \overset{\sim}{m}}^{\frac{1}{4}}},{S_{\overset{\sim}{k}} = {{\eta^{\frac{1}{4}}\tau_{m\quad \overset{\sim}{m}}^{- \frac{1}{4}}} = \tau_{k\overset{\sim}{k}}^{- \frac{1}{4}}}}} & (35)\end{matrix}$

[0123] A point in the parametric space MK{tilde over (K)}{tilde over(M)} is thus completely defined by η and any of these two times. Theevolution of the state point can be conceptualized as moving along atrajectory perpendicular to the storage- or the leak-off-edge.

[0124] In summary, the MK-edge corresponds to the origin of time, andthe {tilde over (M)}{tilde over (K)}-edge to the end of time (except inimpermeable rocks). Thus, given all the problem parameters whichcompletely define the number η, the system evolves with time (e.g., timeτ_(mk)) along a η-trajectory, starting from the MK-edge (C_(m)=C_(k)=0)and ending at the {tilde over (M)}{tilde over (K)}-edge(S_({tilde over (k)})=S_({tilde over (m)})=0). If η=0, the fluid isinviscid (μ′=0) and the system then evolves along the toughness-edgefrom K to {tilde over (K)}. If η=∞, then K′=0 the system evolves alongthe viscosity-edge from M to {tilde over (M)}; The dependence of thescaled solution F can thus be expressed in the form F(ξ;τ;η), where τ isone of the dimensionless time, irrespective of the adopted scaling.

[0125] II. Embodiments Utilizing a Second Parametric Space

[0126] A. Radial Fractures

[0127] Determining the solution of the problem of a radial hydraulicfracture propagating in a permeable rock consists of finding theaperture w of the fracture, and the net pressure p (the differencebetween the fluid pressure p_(f) and the far-field stress σ_(o)) as afunction of both the radial coordinate r and time t, as well as theevolution of the fracture radius R(t). The functions R(t), w(r,t), andp(r,t) depend on the injection rate Q_(o) and on the four materialparameters E′, μ′, K′, and C′ respectively defined as $\begin{matrix}{E^{\prime} = {{\frac{E}{1 - v^{2}}\quad \mu^{\prime}} = {{12\mu \quad K^{\prime}} = {{4\left( \frac{2}{\pi} \right)^{\frac{1}{2}}K_{I\quad c}\quad C^{\prime}} = {2C_{l}}}}}} & (36)\end{matrix}$

[0128] The three functions R(t), w(r,t), and p(r,t) are determined bysolving a set of equations which can be summarized as follows.

[0129] Elasticity Equation $\begin{matrix}{w = {\frac{R}{E^{\prime}}{\int_{0}^{1}{{G\left( {{r/R},s} \right)}{p\left( {{s\quad R},t} \right)}s\quad {s}}}}} & (37)\end{matrix}$

[0130] here G is a known elastic kernel. This singular integral equationexpresses the non-local dependence of the fracture width w on the netpressure p.

[0131] Lubrication Equation $\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{1}{r}\frac{\partial}{\partial r}\left( {r\quad w^{3}\frac{\partial p}{\partial r}} \right)}} & (38)\end{matrix}$

[0132] This non-linear differential equation governs the flow of viscousincompressible fluid inside the fracture. The function g(r,t) denotesthe rate of fluid leak-off, which evolves according to Carter's law$\begin{matrix}{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(r)}}}} & (39)\end{matrix}$

[0133] where t_(o)(r) is the exposure time of point r (i.e., the time atwhich the fracture front was at a distance r from the injection point).

[0134] Global Volume Balance $\begin{matrix}{{Q_{o}t} = {{2\pi {\int_{0}^{R}{w\quad r\quad {r}}}} + {\int_{0}^{t}{r{\int_{0}^{R{(\tau)}}{{g\left( {r,\tau} \right)}\quad {r}\quad {\tau}}}}}}} & (40)\end{matrix}$

[0135] This equation expresses that the total volume of fluid pumped isequal to the sum of the fracture volume and the volume of fluid lost inthe rock surrounding the fracture.

[0136] Propagation Criterion $\begin{matrix}{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{R - r}}},{{1 - \frac{r}{R}}1}} & (41)\end{matrix}$

[0137] Within the framework of linear elastic fracture mechanics, thisequation embodies fact that the fracture is always propagating and thatenergy is dissipated continuously in the creation of new surfaces inrock (at a constant rate per unit surface)

[0138] Tip Conditions

w=0, ${{w^{3}\frac{\partial p}{\partial r}} = 0},$

 r=R  (42)

[0139] The tip of the propagating fracture corresponds to a zero widthand to a zero fluid flow rate condition.

[0140] 1. Scalings

[0141] The general solution of this problem (which includesunderstanding the dependence of the solution on all the problemparameters) can be considerably simplified through the application ofscaling laws. Scaling of this problem hinges on defining thedimensionless crack opening Ω, net pressure Π, and fracture radius γ as

w=εLΩ(ρ;P ₁ ,P ₂), p=εE′Π(ρ;P ₁ ,P ₂), R=γ(P ₁ ,P ₂)L  (43)

[0142] These definitions introduce the scaled coordinate ρ=r/R(t)(0≦ρ≦1), a small number ε(t), a length scale L(t) of the same order ofmagnitude as the fracture length R(t), and two dimensionless evolutionparameters P₁(t) and P₂(t), which depend monotonically on t. As is shownbelow, three different scalings (“viscosity”, “toughness,” and“leak-off”) can be defined, which yield power law dependence of L, ε,P₁, and P₂ on time t; i.e. L˜t^(α), ε˜t^(δ), P₁˜t^(β) ^(₁) , P₂˜t^(β)^(₂) ,. The form of the scaling (43) can be motivated from elementaryelasticity considerations, by noting that the average aperture scaled bythe fracture radius is of the same order as the average net pressurescaled by the elastic modulus.

[0143] The main equations are transformed as follows, under the scaling(43).

[0144] Elasticity Equation $\begin{matrix}{\Omega = {\gamma {\int_{0}^{l}{{G\left( {\rho,s} \right)}{\Pi \left( {{s\text{;}P_{1}},P_{2}} \right)}s{s}}}}} & (44)\end{matrix}$

[0145] Lubrication Equation: $\begin{matrix}{{{\left( {\frac{\overset{.}{ɛ}\quad t}{ɛ} + \frac{\overset{.}{L}\quad t}{L}} \right)\Omega} - {\frac{\overset{.}{L}\quad t}{L}\rho \frac{\partial\Omega}{\partial\rho}} + {{\overset{.}{P}}_{1}{t\left( {\frac{\partial\Omega}{\partial P_{1}} - {\frac{\rho}{\gamma}\frac{\partial\gamma}{\partial P_{1}}\frac{\partial\Omega}{\partial\rho}}} \right)}} + {{\overset{.}{P}}_{2}{t\left( {\frac{\partial\Omega}{\partial P_{2}} - {\frac{\rho}{\gamma}\frac{\partial\gamma}{\partial P_{2}}\frac{\partial\Omega}{\partial\rho}}} \right)}} + {G_{c}\Gamma}} = {\frac{1}{G_{m}\gamma^{2}\rho}\frac{\partial}{\partial\rho}\left( {{\rho\Omega}^{3}\frac{\partial\Pi}{\partial\rho}} \right)}} & (45)\end{matrix}$

[0146] where the leak-off function Γ(ρ;P₁,P₂) is defined as Γ(ρ;P₁,P₂)${{\Gamma \left( {{\rho;P_{1}},P_{2}} \right)} = \frac{1}{\sqrt{1 - {t_{o}/t}}}},$

[0147] t>t_(o)

[0148] Global Mass Balance $\begin{matrix}{{{2\quad {\pi\gamma}^{2}{\int_{0}^{1}{{\Omega\rho}\quad {\rho}}}} + {2\quad \pi \quad G_{c}{\int_{0}^{1}{u^{{2\quad \alpha} - {1/2}}{\gamma^{2}\left( {{u^{\beta_{1}}P_{1}},{u^{\beta_{2}}P_{2}}} \right)}{I\quad\left( {{u^{\beta_{1}}P_{1}},{u^{\beta_{2}}P_{2}}} \right)}{u}}}}} = G_{v}} & (46)\end{matrix}$

[0149] where I is given by I(X₁, X₂) = ∫₀¹Γ(ρ; X₁, X₂)ρ  ρ

[0150] Propagation Criterion

Ω=G _(k)γ^(1/2)(1−ρ)^(1/2)1−ρ<<1  (47)

[0151] Four dimensionless groups G_(v), G_(m), G_(k), G_(c) appear inthese equations: $\begin{matrix}{{G_{v} = \frac{Q_{o}t}{ɛ\quad L^{3}}},{G_{m} = \frac{\mu^{\prime}}{ɛ^{3}E^{\prime}t}},{G_{k} = \frac{K^{\prime}}{ɛ\quad E^{\prime}L^{1/2}}},{G_{c} = \frac{C^{\prime}t^{1/2}}{ɛ\quad L}}} & (48)\end{matrix}$

[0152] While the group G_(v) is associated with the volume of fluidpumped, G_(m), G_(k), and G_(c) can be interpreted as dimensionlessviscosity, toughness, and leak-off coefficients, respectively. Threedifferent scalings can be identified, with each scaling leading to adifferent definition of the set ε, L, P₁, and P₂. Each scaling isobtained by setting G_(v)=1 and one of the other groups to 1 (G_(m) forthe viscosity scaling, G_(k) for the toughness scaling, and G_(c) forthe leak-off scaling), with the two other groups being identified as P₁and P₂. Three scalings denoted as viscosity, toughness, and leak-off canthus be defined depending on whether the group containing μ′ (G_(m)), K′(G_(k)) or C′ (G_(c)) is set to 1. The three scalings are summarized inTable 3. TABLE 3 Small parameter ε, lengthscale L, and parameters P₁ andP₂ for the viscosity, toughness, and leak-off scaling. Scaling ε L P₁ P₂Viscosity $\left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}$

$\left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/9}$

$K_{m} = {K^{\prime}\left( \frac{t^{2}}{\mu^{\prime 5}E^{\prime 13}Q_{o}^{3}} \right)}^{1/18}$

$C_{m} = {C^{\prime}\left( \frac{{E^{\prime}}^{3}t^{7}}{\mu^{\prime 4}Q_{o}^{6}} \right)}^{1/18}$

Toughness$\left( \frac{K^{\prime 6}}{E^{\prime 6}Q_{o}t} \right)^{1/5}$

$\left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/5}$

$M_{k} = {\mu^{\prime}\left( \frac{Q_{o}^{3}E^{\prime 13}}{K^{\prime 18}t^{2}} \right)}^{1/5}$

$C_{k} = {C^{\prime}\left( \frac{E^{\prime 8}t^{3}}{K^{\prime 8}Q_{o}^{2}} \right)}^{1/10}$

Leak-off $\left( \frac{C^{\prime 6}t}{Q_{o}^{2}} \right)^{1/4}$

$M_{\overset{\sim}{k}} = {\mu^{\prime}\left( \frac{E^{\prime 12}C^{\prime 2}Q_{o}^{2}}{K^{\prime 16}t} \right)}^{1/4}$

$K_{c} = {K^{\prime}\left( \frac{Q_{o}^{2}}{E^{\prime 8}C^{\prime 10}t^{3}} \right)}^{1/8}$

$M_{c} = {\mu^{\prime}\left( \frac{Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}t^{7}} \right)}^{1/4}$

[0153] The evolution of the radial fracture can be conceptualized in theternary phase diagram MKC shown in FIG. 10. First, however, thedimensionless number η and time τ are introduced as $\begin{matrix}{{\eta = \frac{K^{\prime 14}}{E^{\prime 11}\mu^{\prime 3}C^{\prime 4}Q_{o}}},{\tau = \frac{t}{t_{c\quad m}}},{{{with}\quad t_{c\quad m}} = \left( \frac{\mu^{\prime 4}Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}} \right)^{1/7}}} & (49)\end{matrix}$

[0154] As shown in Table 3, the evolution parameters P₁ and P₂ in thethree scalings can be expressed in terms of η and τ only. Both K_(m) andC_(m) are positive power of time τ, while K_(c) and M_(c) are negativepower of τ; furthermore, M_(k)˜τ^(−2/5) and C_(k)˜τ^(3/10). Hence, theviscosity scaling is appropriate for small time, while the leak-offscaling is appropriate for large time. The toughness scaling applies tointermediate time when both M_(k) and C_(k) are o(1). TABLE 4 Dependenceof the parameters P₁ and P₂ on the dimensionless time τ and number η forthe viscosity, toughness, and leak-off scaling. Scaling P₁ P₂ ViscosityK_(m) = η^(1/14)τ^(1/9)

C_(m) = τ^(7/18)

Toughness C_(k) = η^(−2/35)τ^(3/10)

M_(k) = η^(−9/35)τ^(−2/5)

Leak-off M_(c) = τ^(−7/4)

K_(c) = η^(1/14)τ^(−3/8)

[0155] The solution of a hydraulic fracture starts at the M-vertex(K_(m)=0, C_(m)=0) and ends at the C-vertex (M_(c)=0, K_(c)=0); itevolves with time τ, along a trajectory which is controlled only by thenumber η, a function of all the problem parameters (i.e., Q_(o), E′, μ′,K′, and C′). If η=0 (the rock has zero toughness), the evolution from Mto C is done directly along the base MC of the ternary diagram MKC. Withincreasing η (which can be interpreted for example as increasingrelative toughness, the trajectory is pulled towards the K-vertex. Forη=∞, two possibilities exist: either the rock is impermeable (C′=0) andthe system evolves directly from M to K, or the fluid is inviscid andthe system then evolves from K to C.

[0156] At each corner of the MKC diagram, there is only one dissipativemechanism at work; for example, at the M-vertex, energy is onlydissipated in viscous flow of the fracturing fluid since the rock isassumed to be impermeable and to have zero toughness. It is interestingto note that the mathematical solution is characterized by a differenttip singularity at each corner, reflecting the different nature of thedissipative mechanism.

[0157] M-Corner:

Ω˜(1−ρ)^(2/3) Π˜(1−ρ)^(−1/3) for ρ˜1  (50)

[0158] C-Corner:

Ω˜(1−ρ)^(5/8) Π˜(1−ρ)^(−3/8) for ρ˜1  (51)

[0159] K-Corner:

Ω˜(1−ρ)^(1/2) Π˜Const for ρ˜1  (52)

[0160] The transition of the solution in the tip region between twocorners can be analyzed by considering the stationary solution of asemi-infinite hydraulic fracture propagating at constant speed.

[0161] 2. Applications of the Scaling Laws

[0162] The dependence of the scaled solution F={Ω,Π,γ} is thus of theform F(ρ,τ;η), irrespective of the adopted scaling. In other words, thescaled solution is a function of the dimensionless spatial and timecoordinates ρ and τ, which depends only on η, a constant for aparticular problem. Thus the laws of similitude between field andlaboratory experiments simply require that η is preserved and that therange of dimensionless time τ is the same—even for the general case whenneither the fluid viscosity, nor the rock toughness, nor the leak-off offracturing fluid in the reservoir can be neglected.

[0163] Although the solution in any scaling can readily be translatedinto another scaling, each scaling is useful because it is associatedwith a particular process. Furthermore, the solution at a corner of theMKC diagram in the corresponding scaling (i.e., viscosity at M,toughness at K, and leak-off at C) is self-similar. In other words, thescaled solution at these vertices does not depend on time, which impliesthat the corresponding physical solution (width, pressure, fractureradius) evolves with time according to a power law. This property of thesolution at the corners of the MKC diagram is important, in part becausehydraulic fracturing near one comer is completely dominated by theassociated process. For example, in the neighborhood of the M-corner,the fracture propagates in the viscosity-dominated regime; in thisregime, the rock toughness and the leak-off coefficient can beneglected, and the solution in this regime is given for all practicalpurposes by the zero-toughness, zero-leak-off solution at the M-vertex.Findings from work along the MK edge where the rock is impermeablesuggest that the region where only one process is dominant is surprisinglarge. FIG. 11 shows the variation of γ_(mk) (the fracture radius in theviscosity scaling) with the dimensionless toughness K_(m) for animpermeable rock (K_(m)=0 corresponds to the M-vertex, K_(m)=∞ (i.e.,M_(c)=0) to the K-vertex). These results indicate that a hydraulicfracture propagating in an impermeable rock is in theviscosity-dominated regime if K_(m)<K_(mm)≅1, and in thetoughness-dominated regime if K_(m)>K_(mk)≅4.

[0164] Accurate solutions can be obtained for a radial hydraulicfracture propagating in regimes corresponding to the edges MK, KC, andCM of the MKC diagram. These solutions enable one to identify the threeregimes of propagation (viscosity, toughness, and leak-off).

[0165] The range of values of the evolution parameters P₁ and P₂ forwhich the fracture propagates in one of the primary regimes (viscosity,toughness, and leak-off) can be identified. The criteria in terms of thenumbers P₁ and P₂ can be translated in terms of the physical parameters(i.e., the injection rate Q_(o), the fluid viscosity μ, the rocktoughness K_(le), the leak-off coefficient C_(l), and the rock elasticmodulus E′).

[0166] The primary regimes of fracture propagation (corresponding to thevertices of the MKC diagram) are characterized by a simple power lawdependence of the solution on time. Along the edges of the MKC triangle,outside the regions of dominance of the corners, the evolution of thesolution can readily be tabulated.

[0167] In some embodiments of the present invention, the tabulatedsolutions are used for quick design of hydraulic fracturing treatments.In other embodiments, the tabulated solutions are used to interpretreal-time measurements during fracturing, such as down-hole pressure.

[0168] The derived solutions can be considered as exact within theframework of assumptions, since they can be evaluated to practically anydesired degree of accuracy. These solutions are therefore usefulbenchmarks to test numerical simulators currently under development.

[0169] 3. Derivation of Solutions Along Edges of the TriangularParametric Space

[0170] Derivation of the solution along the edges of the triangle MKCand at the C-vertex are now described. The identification of thedifferent regimes of fracture propagation are also described.

[0171] a. CK-Solution

[0172] Along the CK-edge of the MKC triangle, the influence of theviscosity is neglected and the solution depends only on one parameter(either K_(c), the dimensionless toughness in the leak-off scaling, orthe dimensionless leak-off coefficient C_(k) in the toughness scalingC_(k)). In one embodiment, the solution is constructed starting from theimpermeable case (K-vertex) and it is evolved with increasing C_(k)towards the C-vertex.

[0173] Since the fluid is taken to be inviscid along the CK-edge, thepressure distribution along the fracture is uniform and thecorresponding opening is directly deduced by integration of theelasticity equation (44) $\begin{matrix}{{\Pi_{k\quad c} = {\Pi_{k\quad c}\left( C_{k} \right)}},{\Omega_{k\quad c} = {\frac{8}{\pi}\gamma_{k\quad c}\Pi_{k\quad c}\sqrt{1 - \rho^{2}}}}} & (53)\end{matrix}$

[0174] Combining (53) with the propagation criterion (47) yields$\begin{matrix}{{\Pi_{k\quad c} = {\frac{\pi}{8\sqrt{2}}\gamma_{k\quad c}^{{- 1}/2}}},{\Omega_{k\quad c} = {\left( \frac{\gamma_{k\quad c}}{2} \right)^{1/2}\sqrt{1 - \rho^{2}}}}} & (54)\end{matrix}$

[0175] The radius γ_(kc) is determined as a function of C_(k). Anequation for γ_(kc) can be deduced from the global balance of mass$\begin{matrix}{{\frac{\gamma_{k\quad c}}{C_{k}} = {{\frac{4}{3\quad C_{k}}{\frac{\gamma_{k}^{5/2}}{\gamma_{k\quad c}^{3/2}}\left\lbrack {1 - \left( \frac{\gamma_{k\quad c}}{\gamma_{k\quad}} \right)^{5/2}} \right\rbrack}} - {\frac{8\quad \pi}{3}\frac{\gamma_{k\quad c}^{1/2}}{\gamma_{k}^{5/2}}{I\left( C_{k} \right)}}}}{where}} & (55) \\{{{\gamma_{k} \equiv {\gamma_{k\quad c}(0)}} = \left( \frac{3}{\pi \sqrt{2}} \right)^{5/2}},{{{and}\quad {I(X)}} = {\int_{0}^{1}{\frac{1}{\sqrt{1 - {\tau_{o}\left( {\rho,X} \right)}}}\quad \rho {\rho}}}}} & (56)\end{matrix}$

[0176] with τ_(o)=t_(o)(r)/t denoting the scaled exposure time of pointr. The function τ_(o)(ρ,X) can be found by inverting $\begin{matrix}{\rho = {\tau_{o}^{2/5}\frac{\gamma_{k\quad c}\left( {\tau_{o}^{3/10}X} \right)}{\gamma_{k\quad c}(X)}}} & (57)\end{matrix}$

[0177] which is deduced from the definition of ρ by taking into accountthe power law dependence of L_(k) and C_(k) on time.

[0178] Since τ_(o)(1,X)=1, the integral I(X) defined in (56) is singularat ρ=1. This singularity is weak, and its strength is known at X=0 andX=∞: X=0(τ_(o)=ρ^(5/2)) and at X=∞(τ_(o)=τ⁴). From a computational pointof view, the integral can be calculated along the time axis with respectto τ_(o) $\begin{matrix}{{I(X)} = {\frac{1}{\gamma_{k\quad c}(X)}{\int_{0}^{1}{{\frac{1}{{\tau_{o}^{3/5}\left( {1 - \tau_{o}} \right)}^{1/2}}\left\lbrack {{\frac{2}{5}{\gamma_{k\quad c}\left( {\tau_{o}^{3/10}X} \right)}} + {\frac{3}{10}\tau_{o}^{3/10}X\quad {\gamma_{k\quad c}^{\prime}\left( {\tau_{o}^{3/10}X} \right)}}} \right\rbrack}{\tau_{o}}}}}} & (58)\end{matrix}$

[0179] In some embodiments of the present invention, the solution can beobtained by solving the non-linear ordinary differential equation (55),using an implicit iterative algorithm.

[0180] b. MK-Solution

[0181] The MK-solution corresponds to regimes of fracture propagation inimpermeable rocks. One difficulty in obtaining this solution lies inhandling the changing nature of the tip behavior between the M- and theK-vertex. The tip asymptote is given by the classical square rootsingularity of linear elastic fracture mechanics (LEFM) wheneverK_(m)≠0. However, near the M-vertex (small K_(m)), the LEFM behavior isconfined to a small boundary layer, which does not influence thepropagation of the fracture. In this viscosity-dominated regime, thesingularity (50) develops as an intermediate asymptote.

[0182] The solution can be searched for in the form of a finite seriesof known base functions $\begin{matrix}{\Pi_{k\quad m} = {{\Pi_{o}^{*}\left( {\rho,M_{k}} \right)} + {\sum\limits_{i = 1}^{n_{\Pi}}\quad {{A_{i}\left( M_{k} \right)}{\Pi_{i}^{*}(\rho)}}} + {{B\left( M_{k} \right)}{\Pi^{**}(\rho)}}}} & (59) \\{{\overset{\_}{\Omega}}_{k\quad m} = {{{\overset{\_}{\Omega}}_{o}^{*}\left( {\rho,M_{k}} \right)} + {\sum\limits_{i = 1}^{n_{\Omega}}\quad {{C_{i}\left( M_{k} \right)}{{\overset{\_}{\Omega}}_{i}^{*}(\rho)}}} + {{B\left( M_{k} \right)}{{\overset{\_}{\Omega}}^{**}(\rho)}}}} & (60)\end{matrix}$

[0183] where the introduction of {overscore (Ω)}_(km)=Ω_(km)/γ_(km)excludes γ_(km) from the elasticity equation (44).

[0184] Since the non-linearity of the problem only arises from thelubrication equation (45), the series expansions (59) and (60) can beused to satisfy the elasticity equation and the boundary conditions atthe tip and at the inlet. In the proposed decomposition, the last terms{Π**,{overscore (Ω)}**} are chosen such that the logarithmic pressuresingularity near the inlet is satisfied. The corresponding opening isintegrated by substituting this pressure function into (44). The firstterms in the series {Π_(o)*,{overscore (Ω)}_(o)*} are constructed toexactly satisfy the propagation equation and to account for thelogarithmic pressure asymptote near the tip (which results fromsubstituting the opening square root asymptote into the lubricationequation). It is also required that {Π_(o)*,{overscore (Ω)}_(o)*}exactly satisfy the elasticity equation (44). The regular part of thesolution is represented by series of base functions {Π_(i)*,{overscore(Ω)}_(i)*} The choice of these functions is not unique; however, itseems consistent to require that {overscore (Ω)}_(i)*˜(1−ρ)^(1/2+i) forρ˜1. (The square root opening asymptote appears only in the first term,if one imposes that the function Π_(i)*; does not contribute to thestress intensity factor.) A convenient choice of these base functionsare Jacobi polynomials with the appropriate weights.

[0185] Any pair {Π_(i)*,{overscore (Ω)}_(i)*} does not satisfy theelasticity equation (44). Instead, the coefficients A_(i) and C_(i) arerelated by the elasticity equation through the matrix L_(ij) (which isindependent of M_(k) or K_(m)). $\begin{matrix}{C_{i}^{({n_{\Omega},n_{\Pi}})} = {\sum\limits_{j = 1}^{n_{\Pi}}{L_{i\quad j}A_{j}^{({n_{\Omega},n_{\Pi}})}}}} & (61)\end{matrix}$

[0186] The problem is reduced to finding n_(Π)+1 unknown coefficientsA_(i) and B, by solving the lubrication equation (45), which simplifieshere to $\begin{matrix}\begin{matrix}{{\rho {\overset{\_}{\Omega}}_{k\quad m}^{3}\frac{\partial\Pi_{k\quad m}}{\partial\rho}} + {M_{k}{\int_{\rho}^{1}{{\overset{\_}{\Omega}}_{k\quad m}s\quad {s}}}} + {\frac{4}{5}M_{k}\rho^{2}{\overset{\_}{\Omega}}_{k\quad m}} -} \\{{{- \frac{2}{5}}{M_{k}^{2}\left\lbrack {{\int_{\rho}^{1}{\frac{\partial{\overset{\_}{\Omega}}_{k\quad m}}{\partial M_{k}}s\quad {s}}} + {\frac{2}{\gamma_{k\quad m}}\frac{\gamma_{k\quad m}}{M_{k}}\left( {{\int_{\rho}^{1}{{\overset{\_}{\Omega}}_{k\quad m}s\quad {s}}} + {\rho^{2}{\overset{\_}{\Omega}}_{k\quad m}}} \right)}} \right\rbrack}} = 0}\end{matrix} & (62)\end{matrix}$

[0187] where γ_(km)=(2π∫₉₂ {overscore (Ω)}_(km)ρdρ)^(−1/3)

[0188] In some embodiments of the present invention, the lubricationequation is solved by an implicit iterative procedure. For example, thesolution at the current iteration can be found by a least squaresmethod.

[0189] c. CM-Solution

[0190] In some embodiments, the solution along the CM-edge of the MKCtriangle is found using the series expansion technique described abovewith reference to the MK-solution. In other embodiments, a numericalsolution is used based on the following algorithm.

[0191] The displacement discontinuity method is used to solve theelasticity equation (44). This method yields a linear system ofequations between aperture and net pressure at nodes along the fracture.The coefficients (which can be evaluated analytically) need to becalculated only once as they do not depend on C_(m). The lubricationequation (45) is solved by a finite difference scheme (either explicitor implicit). The fracture radius γ_(mc) is found from the global massbalance. Here, the numerical difficulty is to calculate the amount offluid lost due to the leak-off.

[0192] The propagation is governed by the asymptotic behavior of thesolution at the fracture tip. The tip asymptote can be used to establisha relationship between the opening at the computational node next to thetip and the tip velocity. However, this relationship evolves as C_(m)increases from 0 to ∞ (i.e., when moving from the M- to the C-vertex);it is obtained through a mapping of the autonomous solution of asemi-infinite hydraulic fracture propagating at constant speed in apermeable rock.

[0193] d. Solution Near the C-Vertex

[0194] The limit solution at the C-vertex, where both the viscosity andthe toughness are neglected, is degenerated as all the fluid injectedinto the fracture has leaked into the rock. Thus the opening and the netpressure of the fracture is zero, while its radius is finite. In someembodiments of the present invention, the solution near the C-vertex isused for testing the numerical solutions along the CK and CM sides ofthe parametric triangle. The limitation of those solutions comes fromthe choice of the scaling. In order to approach the C-vertex, thecorresponding parameter (C_(k) or C_(m)) must grow indefinitely.Practically, these solutions are calculated up to some finite values ofthe parameters, for which they can be connected with asymptoticsolutions near the C-vertex along CM and CK sides. These asymptoticsolutions can be constructed as follows.

[0195] Consider first the CM-solutionF_(cm)={Ω_(cm)(ρ,M_(c)),Π_(cm)(ρ,M_(c)),γ_(cm)(M_(c))} near theC-vertex. It can be asymptotically approximated as

γ_(cm)=γ_(c) o(M _(c)), Ω_(cm) =M _(c) ^(α)γ_(c){overscore(Ω)}_(cm)(ρ)+o(M _(c) ^(α)), Π_(cm) =M _(c) ^(α){overscore(Π)}_(cm)(ρ)+o(M _(c) ^(α))  (63)

[0196] where γ_(c) denotes the finite fracture radius (in the leak-offscaling) at the C-vertex. The exponent α=¼ is determined by substitutingthese expansions into the lubrication equation (45), which then reducesto $\begin{matrix}{\frac{\gamma_{c}}{\sqrt{1 - \rho^{4}}} = {\frac{1}{\rho}\frac{\quad}{\rho}\left( {\rho {\overset{\_}{\Omega}}_{c\quad m}^{3}\frac{{\overset{\_}{\Pi}}_{c\quad m}}{\rho}} \right)}} & (64)\end{matrix}$

[0197] The asymptotic solution {overscore (F)}_(cm)={{overscore(Ω)}_(cm)(ρ),{overscore (Π)}_(cm)(ρ)} near the C-vertex is found bysolving (64) along with the elasticity equation (44). This can be doneusing the series expansion technique described above. This problem issimilar to the problem at the M-vertex (fracture propagating in animpermeable solid with zero toughness), but with a different tipasymptote. Thus a set of base functions different from the one used forthe M-corner are introduced.

[0198] The CK-solutionF_(ck)={Ω_(ck)(ρ,K_(c)),Π_(ck)(ρ,K_(c)),γ_(ck)(K_(c))} near the C-vertexcan also be sought in the form of an asymptotic expansion${\gamma_{c\quad k} = {\gamma_{c} + {o\left( K_{c} \right)}}},{\Omega_{ck} = {{K_{c}^{\beta}\gamma_{c}{{\overset{\_}{\Omega}}_{ck}(\rho)}} + {o\left( K_{c}^{\beta} \right)}}},{\Pi_{ck} = {{K_{c}^{\beta}{{\overset{\_}{\Pi}}_{c\quad k}(\rho)}} + {o\left( K_{c}^{\beta} \right)}}}$

[0199] where β=1 is determined from the propagation condition (11). Thissolution is trivial, however, since the pressure is uniform; hence$\begin{matrix}{{{\overset{\_}{\Pi}}_{c\quad k} = {\frac{\pi}{8}\left( {2\quad \gamma_{c}} \right)^{{- 1}/2}}},{{\overset{\_}{\Omega}}_{ck} = {\left( {2\gamma_{c}} \right)^{{- 1}/2}\sqrt{1 - \rho^{2}}}}} & (66)\end{matrix}$

[0200] e. Regimes of Fracture Propagation

[0201] The regimes of fracture propagation can readily be identifiedonce the solutions at the vertices and along the edges of the MKCtriangle have been tabulated using the algorithms and methods ofsolutions described above. Recall that for the parametric space underconsideration, there are three primary regimes of propagation(viscosity, toughness, and leak-off) associated with the vertices of theMKC triangle and that in a certain neighborhood of a corner, thecorresponding process is dominant, see Table 5. For example, fracturepropagation is in the viscosity-dominated regime if K_(m)<K_(mm) andC_(m)<C_(mm); in this region, the solution can be approximated for allpractical purposes by the zero-toughness, zero-leak-off solution at theM-corner (K_(m)=0, C_(m)=0) TABLE 5 Range of the parameters P₁ and P₂for which a primary process is dominant. Dominant Process Range on P₁Range on P₂ Viscosity K_(m) < K_(mm) (M_(k) > M_(km)) C_(m) < C_(mm)(M_(c) > M_(cm)) Toughness C_(k) < C_(kk) (K_(c) > K_(ck)) M_(k) <M_(kk) (K_(m) > K_(mk)) Leak-off M_(c) < M_(cc) (C_(m) > C_(mc)) K_(c) <K_(cc) (C_(k) < C_(kc))

[0202] Identification of the threshold values of the evolutionparameters (for example, K_(mm) and C_(mm) for the viscosity-dominatedregime) can be accomplished by comparing the fracture radius with itsreference value at a corner. The corner process is considered asdominant, if the fracture radius is within 1% of its value at thecorner. For example, K_(mm) and C_(mm) are deduced from the followingconditions

|γ_(mk)(K _(mm))−γ_(m)|/γ_(m)=1%, |γ_(mk)(C _(mm))−γ_(m)|/γ_(m)=1%  (67)

[0203] B. Plane Strain (KGD) Fractures

[0204] 1. Governing Equations and Boundary Conditions

[0205] Elasticity $\begin{matrix}{w = {\frac{l}{E^{\prime}}{\int_{0}^{1}{{G\left( {{x/l},s} \right)}{\Pi \left( {{s\quad l},t} \right)}\quad {s}}}}} & (68)\end{matrix}$

[0206] Lubrication $\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{\partial\quad}{\partial x}\left( {w^{3}\frac{\partial p}{\partial x}} \right)}} & (69)\end{matrix}$

[0207] obtained by eliminating the radial flow rate q(x,t) between thefluid mass balance $\begin{matrix}{{\frac{\partial w}{\partial t} + \frac{\partial q}{\partial x} + g} = 0} & (70)\end{matrix}$

[0208] and Poiseuille law $\begin{matrix}{q = {{- \frac{w^{3}}{\mu^{\prime}}}\frac{\partial p}{\partial x}}} & (71)\end{matrix}$

[0209] Leak-Off $\begin{matrix}{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(x)}}}} & (72)\end{matrix}$

[0210] where t_(o)(x) is the exposure time of point x

[0211] Global Volume Balance $\begin{matrix}{{Q_{o}t} = {{2{\int_{0}^{l{(t)}}{w\quad {x}}}} + {2{\int_{0}^{l}{\int_{0}^{l{(\tau)}}\quad {{g\left( {x,\tau} \right)}{x}{\tau}}}}}}} & (73)\end{matrix}$

[0212] Propagation Criterion $\begin{matrix}{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{l - x}}},{{1 - \frac{x}{l}}1}} & (74)\end{matrix}$

[0213] Tip Conditions

w=0, $\begin{matrix}{{{w^{3}\frac{\partial p}{\partial x}} = 0},\quad {x = l}} & (75)\end{matrix}$

[0214] 2. Scaling

[0215] Similarly to the radial fracture, we define the dimensionlesscrack opening Ω, net pressure Π, and fracture length γ as

w(x,t)=ε(t)L(t)Ω(ξ;P ₁ P ₂)  (76)

p(x,t)=ε(t)E′Π(ξ;P ₁ P ₂)  (77)

l(t)=γ(ξ;P ₁ P ₂)L(t)  (78)

[0216] These definitions introduce a scaled coordinate ξ=x/l(t) (0≦ξ≦1),a small number ε(t), a length scale L(t) of the same order of magnitudeas the fracture length l(t), and two dimensionless parameters P₁(t),P₂(t) which depend monotonically on t. The form of the scaling (76)-(80)can be motivated from elementary elasticity considerations, by notingthat the average aperture scaled by the fracture radius is of the sameorder as the average net pressure scaled by the elastic modulus.Explicit forms of the parameters ε(t), L(t), P₁(t), and P₂ (t) are givenbelow for the viscosity, toughness, and leak-off scalings.

[0217] The main equations are transformed as follows, under the scaling(76)-(80).

[0218] Elasticity Equation

Ω=γ∫₀ ¹ G(ξ,s)Π(s;P ₁ ,P ₂)ds  (79)

[0219] Lubrication Equation.

[0220] The left-hand side ∂w/∂t of the lubrication equation (69) can nowbe written as $\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {{\left( {{\overset{.}{ɛ}\quad L} + {ɛ\quad \overset{.}{L}}} \right)\Omega} - {ɛ\quad \overset{.}{L}\xi \frac{\partial\Omega}{\partial\xi}} + {ɛ\quad L\quad {{\overset{.}{P}}_{1}\left( {\frac{\partial\Omega}{\partial P_{1}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{1}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {ɛ\quad L\quad {{\overset{.}{P}}_{2}\left( {\frac{\partial\Omega}{\partial P_{2}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{2}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {C^{\prime}t^{- \frac{1}{2}}{\Gamma \left( {{\xi \text{;}P_{1}},P_{2}} \right)}}}} & (80)\end{matrix}$

[0221] while the right hand side is transformed into $\begin{matrix}{{\frac{1}{\mu^{\prime}}\frac{\partial p}{\partial x}\left( {w^{3}\frac{\partial p}{\partial x}} \right)} = {\frac{ɛ^{4}E^{\prime}L}{\mu^{\prime}\gamma^{2}}\frac{\partial}{\partial\xi}\left( {\Omega^{3}\frac{\partial\Pi}{\partial\xi}} \right)}} & (81)\end{matrix}$

[0222] The leak-off function Γ(ξ;P₁,P₂), which is defined as$\begin{matrix}{{{\Gamma \left( {{\xi \text{;}P_{1}},P_{2}} \right)} = \frac{1}{\sqrt{1 - {t_{o}/t}}}},{t > t_{o}}} & (82)\end{matrix}$

[0223] can be computed as part of the solution, once the parameters P₁,P₂ have been identified. After multiplying both sides by t/εR, we obtaina new form of the lubrication equation $\begin{matrix}{{{\left( {\frac{\overset{.}{ɛ}\quad t}{ɛ} + \frac{\overset{.}{L}\quad t}{L}} \right)\Omega} - {\frac{\overset{.}{L}\quad t}{L}\xi \frac{\partial\Omega}{\partial\xi}} + {{\overset{.}{P}}_{1}{t\left( {\frac{\partial\Omega}{\partial P_{1}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{1}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {{\overset{.}{P}}_{2}{t\left( {\frac{\partial\Omega}{\partial P_{2}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{2}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {\frac{C^{\prime}t^{\frac{1}{2}}}{ɛ\quad L}\Gamma}} = {\frac{ɛ^{3}E^{\prime}t}{\mu^{\prime}}\frac{1}{\gamma^{2}}\frac{\partial}{\partial\xi}\left( {\Omega^{3}\frac{\partial\Pi}{\partial\xi}} \right)}} & (83)\end{matrix}$

[0224] Global Mass Balance $\begin{matrix}{{{2\gamma {\int_{0}^{l}{\Omega \quad {\rho}}}} + {2\frac{C^{\prime}t^{\frac{1}{2}}}{ɛ\quad L}{\int_{0}^{l}{u^{\alpha - {1/2}}{\gamma \left( {{\tau_{o}^{\beta_{1}}P_{1}},{\tau_{o}^{\beta_{2}}P_{2}}} \right)}{I\left( {{\tau_{o}^{\beta_{1}}P_{1}},{\tau_{o}^{\beta_{2}}P_{2}}} \right)}{u}}}}} = \frac{Q_{o}t}{ɛ\quad L^{2}}} & (84)\end{matrix}$

[0225] where I is given by I(X₁,X₂)=∫₀ ¹Γ(ξ;X₁,X₂)dξ

[0226] Propagation Criterion $\begin{matrix}{\Omega = {{{\frac{K^{\prime}}{ɛ\quad E^{\prime}\quad L^{\frac{1}{2}}}{\gamma^{\frac{1}{2}}\left( {1 - \xi} \right)}^{\frac{1}{2}}1} - \xi}1}} & (85)\end{matrix}$

[0227] These equations show that there are 4 dimensionless groups:G_(v), G_(m), G_(k), G_(c), (only G_(v) differs from the radial case, inview of the different dimension of Q_(o)) $\begin{matrix}{{G_{v} = \frac{Q_{o}t}{ɛ\quad L^{2}}},{G_{m} = \frac{\mu^{\prime}}{ɛ^{3}E^{\prime}t}},{G_{k} = \frac{K^{\prime}}{ɛ\quad E^{\prime}L^{1/2}}},{G_{c} = \frac{C^{\prime}t^{1/2}}{ɛ\quad L}}} & (86)\end{matrix}$

[0228] a. Viscosity Scaling.

[0229] The small parameter ε_(m) and the lengthscale L_(m) aredetermined by setting G_(v)=1 and G_(m)=1. Hence, $\begin{matrix}{{ɛ_{m} = \left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}},{L_{m} = \left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/6}}} & (87)\end{matrix}$

[0230] The two parameters P₁=G_(k) and P₂=G_(c) are identified as K_(m)and C_(m), a dimensionless toughness and a dimensionless leak-offcoefficient, respectively $\begin{matrix}{{K_{m} = {K^{\prime}\left( \frac{1}{E^{\prime 3}\mu^{\prime}Q_{o}} \right)}^{1/4}},{C_{m} = {C^{\prime}\left( \frac{E^{\prime}t}{\mu^{\prime}Q_{o}^{3}} \right)}^{1/6}}} & (88)\end{matrix}$

[0231] b. Toughness Scaling.

[0232] Now, ε_(k) and L_(k) are determined from G_(v)=1 and G_(k)=1.Hence, $\begin{matrix}{{ɛ_{k} = \left( \frac{K^{\prime 4}}{E^{\prime 4}Q_{o}t} \right)^{1/3}},{L_{k} = \left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/3}}} & (89)\end{matrix}$

[0233] The two parameters P₁=G_(m) and P₂=G_(c) correspond to M_(k) andC_(k), a dimensionless viscosity and a dimensionless leak-offcoefficient, respectively $\begin{matrix}{{M_{k} = {\mu^{\prime}\left( \frac{E^{\prime 3}Q_{o}}{K^{\prime 4}} \right)}},{C_{k} = {C^{\prime}\left( \frac{E^{\prime 4}t}{K^{\prime 4}Q_{o}^{2}} \right)}^{1/6}}} & (90)\end{matrix}$

[0234] c. Leak-Off Scaling.

[0235] Finally, the leak-off scaling corresponds to G_(v)=1 and G_(c)=1.Hence, $\begin{matrix}\begin{matrix}{ɛ_{c} = \frac{C^{\prime 2}}{Q_{o}}} & {L_{c} = \left( \frac{Q_{o}^{2}t}{C^{\prime 2}} \right)^{1/2}}\end{matrix} & (91)\end{matrix}$

[0236] and the two parameters P₁=G_(k) and P₂=G_(m) are now identifiedas K_(c) and M_(c), a dimensionless viscosity and a dimensionlesstoughness, respectively $\begin{matrix}{{K_{c} = {K^{\prime}\left( \frac{Q_{o}^{2}}{E^{\prime 4}C^{\prime 6}t} \right)}^{1/4}},{M_{c} = {\mu^{\prime}\left( \frac{Q_{o}^{3}}{E^{\prime}C^{\prime 6}t} \right)}}} & (92)\end{matrix}$

[0237] We note that both C_(k), C_(m) are positive power of time t whileK_(c), M_(c) are negative power of t. Hence, the leak-off scaling isappropriate for large time, and either the viscosity scaling or thetoughness scaling is appropriate for small time. As discussed below, thesolution starts from a point on the MK-side of a ternary parameter space(C_(k)=0, C_(m)=0) and tends asymptotically towards the C-point(M_(c)=0, K_(c)=0), following a straight trajectory which is controlledby a certain number η, a function of all the problem parameters exceptC′ (i.e., Q_(o), E′, μ′, K′)

[0238] 3. Time Scales

[0239] It is of interest to express the small parameters ε's, the lengthscales L's, and the dimensionless parameters M's, K's, and C's in termsof time scales. Two time scales t_(m), t_(k) are naturally defined as$\begin{matrix}{{t_{m} = \frac{\mu^{\prime}}{E^{\prime}}},{t_{k} = \frac{K^{\prime 4}}{E^{\prime 4}Q_{o}}}} & (93)\end{matrix}$

[0240] Note that unlike the radial fracture, it is not possible todefine a characteristic time t_(c), since Q_(o) has the dimensionsquared of C′. Hence, $\begin{matrix}{{ɛ_{m} = \left( \frac{t_{m}}{t} \right)^{1/3}},{ɛ_{k} = \left( \frac{t_{k}}{t} \right)^{1/3}}} & (94) \\{{L_{m} = {\left( \frac{t}{t_{m}} \right)^{2/3}{\overset{\_}{L}}_{m}}},{L_{k} = {\left( \frac{t}{t_{k}} \right)^{2/3}{\overset{\_}{L}}_{k}}}} & (95)\end{matrix}$

[0241] where the {overscore (L)}'s are intrinsic length scales definedas $\begin{matrix}{{{\overset{\_}{L}}_{m} = \left( \frac{\mu^{\prime 3}Q_{o}^{3}}{E^{\prime 3}} \right)^{1/6}},{{\overset{\_}{L}}_{k} = \left( \frac{K^{\prime}}{E^{\prime}} \right)^{2}}} & (96)\end{matrix}$

[0242] Next, consider the dimensionless parameters M's, K's, and C'swhich can be rewritten in terms of the characteristic times t_(cm), andt_(ck) $\begin{matrix}{{C_{m} = \left( \frac{t}{t_{cm}} \right)^{1/6}},\quad {C_{k} = \left( \frac{t}{t_{ck}} \right)}} & (97) \\{{M_{c} = \left( \frac{t_{cm}}{t} \right)},\quad {K_{c} = \left( \frac{t_{ck}}{t} \right)^{1/4}}} & (98) \\{{{{where}\quad t_{cm}} = {{ɛ_{c}^{- 3}t_{m}} = \left( \frac{\mu^{\prime}Q_{o}^{3}}{E^{\prime}C^{\prime 6}} \right)}},\quad {t_{ck} = {{ɛ_{c}^{- 3}t_{k}} = \left( \frac{K^{\prime 4}Q_{o}^{2}}{E^{\prime 4}C^{\prime 6}} \right)}}} & (99)\end{matrix}$

[0243] It is thus convenient to introduce a parameter η related to theratios of characteristic times, which is defined as $\begin{matrix}{\eta = \frac{K^{\prime 4}}{E^{\prime 3}\mu^{\prime}Q_{o}}} & (100)\end{matrix}$

[0244] Indeed, it is easy to show that the various characteristic timeratios can be expressed in terms of η $\begin{matrix}{\frac{t_{ck}}{t_{cm}} = \eta} & (101)\end{matrix}$

[0245] Note also that η can be expressed as $\begin{matrix}{\eta = \frac{{\overset{\_}{L}}_{k}^{2}}{{\overset{\_}{L}}_{m}^{2}}} & (102)\end{matrix}$

[0246] Furthermore, if we introduce the dimensionless time τ$\begin{matrix}{\tau = \frac{t}{t_{cm}}} & (103)\end{matrix}$

[0247] (acknowledging at the same time that the choice of t_(cm) toscale the time is arbitrary, as t_(ck) could have been used as well),the parameters M's, K's, and C's can be expressed in terms of τ and η asfollows

K _(m)=η^(1/4) , C _(m)=τ^(1/6) , C _(k)η=^(−1/6)τ^(1/6)  (104)

M _(k)=η⁻¹ , M _(c)=τ⁻¹ , K _(c)=η^(1/4)τ^(−1/4)  (105)

[0248] The dependence of the scaled solution F={Ω,Π,γ} is thus of theform F(ξ,τ;η), irrespective of the adopted scaling (but γ=γ(τ;η)). Inother words, the scaled solution is a function of dimensionless spatialand time coordinate, ξ and τ, which depends on only one parameter, η,which is constant for a particular problem. Thus the laws of similitudebetween field and laboratory experiments simply require that η ispreserved and that the range of dimensionless time τ is the same—evenfor the general case of viscosity, toughness, and leak-off.

[0249] It is again convenient to introduce the ternary diagram MKC shownin FIG. 12. With time τ, the system evolves along a η-trajectory (alongwhich η is a constant), starting from a point on the MK-side and endingat the C-vertex. If η=0 (the rock has zero toughness), the evolutionfrom M to C is done directly along the base BC of the ternary diagramMKC. For η=∞, the fluid is inviscid and the system then evolves from Kto C.

[0250] The KGD fracture differs from the radial fracture by theexistence of only characteristic time rather than two for thepenny-shaped fracture. The characteristic number η for the KGD fractureis independent of the leak-off coefficient C′, which only enters thescaling of time.

[0251] 4. Relationship Between Scalings

[0252] Any scaling can be translated into any of the other two. It canreadily be established that

K _(m) =M _(k) ^(−1/4) , C _(m) =M _(c) ^(−1/6) , C _(k) =K _(c)^(−2/3)  (106)

[0253] and $\begin{matrix}{\frac{ɛ_{m}}{ɛ_{k}} = {M_{k}^{1/3} = K_{m}^{{- 4}/3}}} & (107) \\{\frac{ɛ_{k}}{ɛ_{c}} = {K_{c}^{4/3} = C_{k}^{- 2}}} & (108) \\{\frac{ɛ_{c}}{ɛ_{m}} = {C_{m}^{2} = M_{c}^{{- 1}/3}}} & (109) \\{\frac{L_{m}}{L_{k}} = {M_{k}^{{- 1}/6} = K_{m}^{2/3}}} & (110) \\{\frac{L_{k}}{L_{c}} = {K_{c}^{{- 2}/3} = C_{k}}} & (111) \\{\frac{L_{c}}{L_{m}} = {C_{m}^{- 1} = M_{c}^{1/6}}} & (112)\end{matrix}$

[0254] III. Applications

[0255] Applications of hydraulic fracturing include the recovery of oiland gas from underground reservoirs, underground disposal of liquidtoxic waste, determination of in-situ stresses in rock, and creation ofgeothermal energy reservoirs. The design of hydraulic fracturingtreatments benefits from information that characterize the fracturingfluid, the reservoir rock, and the in-situ state of stress. Some ofthese parameters are easily determined (such as the fluid viscosity),but for others, it is more difficult (such as physical parameterscharacterizing the reservoir rock and in-situ state of stress).

[0256] By utilizing the various embodiments of the present invention,the “difficult” parameters can be assessed from measurements (such asdownhole pressure) collected during a hydraulic fracturing treatment.The various embodiments of the present invention recognize that scaledmathematical solutions of hydraulic fractures with simple geometrydepend on only two numbers that lump time and all the physicalparameters describing the problem. There are many different ways tocharacterize the dependence of the solution on two numbers, as describedin the different sections above, and all of these are within the scopeof the present invention.

[0257] Various parametric spaces have been described, and trajectorieswithin those spaces have also been described. Each trajectory shows apath within the corresponding parametric space that describes theevolution of a particular treatment over time for a given set ofphysical parameter values. That is to say, each trajectory lumps all ofthe physical parameters, except time. Since there exists a uniquesolution at each point in a given parametric space, which needs to becalculated only once and which can be tabulated, the evolution of thefracture can be computed very quickly using these pre-tabulatedsolutions. In some embodiments, pre-tabulated points are very closetogether in the parametric space, and the closest pre-tabulated point ischosen as a solution. In other embodiments, solutions are interpolatedbetween pre-tabulated points.

[0258] The various parametric spaces described above are useful toperform parameter identification, also referred to as “data inversion.”Data inversion involves solving the so-called “forward model” manytimes, where the forward model is the tool to predict the evolution ofthe fracture, given all the problems parameters. Data inversion alsoinvolves comparing predictions from the forward model with measurements,to determine the set of parameters that provide the best match betweenpredicted and measured responses.

[0259] Historically, running forward models has been computationallydemanding, especially given the complexity of the hydraulic fracturingprocess. Utilizing the various embodiments of the present invention,however, the forward model includes pre-tabulated scaled solutions interms of two dimensionless parameters, which only need to be “unsealed”through trivial arithmetic operations. These developments, and others,make possible real-time, or near real-time, data inversion whileperforming a hydraulic fracturing treatment.

[0260] Although the present invention has been described in conjunctionwith certain embodiments, it is to be understood that modifications andvariations may be resorted to without departing from the spirit andscope of the invention as those skilled in the art readily understand.Such modifications and variations are considered to be within the scopeof the invention and the appended claims. For example, the scope of theinvention encompasses the so-called power law fluids (a generalizedviscous fluid characterized by two parameters and which degenerates intoa Newtonian fluid when the power law index is equal to 1). Also forexample, the scope of the invention encompasses the evolution of thehydraulic fracture following “shut-in” (when the injection of fluid isstopped). Hence, various embodiments of the invention contemplateinterpreting data gathered after shut-in.

What is claimed is:
 1. A method comprising: receiving hydraulicfracturing treatment data; and evaluating a forward model to predict theevolution of a fracture, wherein the forward model comprisespre-tabulated scaled solutions in terms of at least one dimensionlessparameter.
 2. The method of claim 1 further comprising unscaling thepre-tabulated solutions to produce a value for at least one physicalparameter.
 3. The method of claim 1 wherein the forward model comprisespre-tabulated solutions in terms of at least two dimensionless evolutionparameters.
 4. The method of claim 3 wherein one of the dimensionlessparameters represents a dimensionless leak-off coefficient.
 5. Themethod of claim 3 wherein the at least two dimensionless evolutionparameters comprise monotonic functions of time.
 6. The method of claim1 wherein the hydraulic fracturing treatment data comprises a pressureof a viscous fluid.
 7. The method of claim 1 wherein the hydraulicfracturing treatment data comprises a fracture dimension.
 8. A methodcomprising: injecting a viscous fluid to fracture reservoir rock;collecting data from measurements made during the injecting; andidentifying values of parameters that characterize the reservoir rockfrom the data, wherein identifying comprises unscaling pre-tabulatedsolutions in terms of at least one dimensionless parameter.
 9. Themethod of claim 8 further comprising identifying a value of a parameterthat characterizes a toughness of the reservoir rock by unscaling apre-tabulated solution.
 10. The method of claim 8 wherein identifyingcomprises identifying a leak-off coefficient.
 11. The method of claim 8wherein identifying comprises identifying a rock permeability.
 12. Amethod of hydraulic fracturing comprising: injecting a viscous fluid;measuring a pressure of the viscous fluid; determining a value of atleast one dimensionless parameter associated with the pressure; anddetermining a value of a physical parameter from the at least onedimensionless parameter.
 13. The method of claim 12 wherein at least oneof the dimensionless parameters represents a toughness of reservoirrock.
 14. The method of claim 12 wherein at least one of thedimensionless parameters represents a fluid storage mechanism.
 15. Themethod of claim 12 wherein at least one of the dimensionless parametersrepresents a dimensionless leak-off coefficient.
 16. The method of claim12 wherein injecting a viscous fluid comprises injecting a viscous fluidat a substantially constant volumetric rate.
 17. A method of designing ahydraulic fracturing treatment comprising: storing pre-tabulatedsolutions that represent problem solution points in a parametric space,wherein the parametric space corresponds to a scaling of the problemparameters; and determining a volumetric rate based on a desiredtrajectory in the parametric space.
 18. The method of claim 17 whereinthe scaling comprises a dimensionless crack opening.
 19. The method ofclaim 17 wherein the scaling comprises a dimensionless net pressure. 20.The method of claim 17 wherein the scaling comprises a dimensionlessfracture radius.